1. Introduction

2. Theoretical Background of Gradient Elasticity

3. Smoothed finite element method: overview

4. Numerical tests

4.1 Cantilever beam

4.2 Plate with a U-notched hole

4.3 Plate with a crack

5. Conclusions

^{} 1. Introduction

Classical continuum mechanics theories are widely used in various engineering applications. Through all-embracing experiments, it has been found that many engineering materials, e.g. metals, composite, concrete, ceramics, polymers, etc., behave differently when they are of different sizes. Such mechanical behaviors by the size effect have been revealed by the micro-structure of the materials (Kolo, 2019). Classical elasticity theory troubles capturing the mechanical behavior of the size effect often observed at a micro-scale level because their assumption is obviously larger than micro-scale. Therefore, this leads to the development of gradient elasticity which possesses the internal length scale (Tenek and Aifantis, 2001; 2002). The internal length scale represents the information on the micro-structure of materials. However, its physical interpretation still remains with no consensus since it depends on material context.

To properly explain it, various studies have been proposed: strain gradient theory (Aifantis, 1992; Ru and Aifantis, 1993), gradient continua (Auffray *et al*., 2015), modified couple stress theory (Mindlin, 1964; Toupin, 1962) and non-local elasticity theory (Eringen, 1983; Pisano *et al*., 2009). These methods have one or more additional constants. Hence, one of the simplest theories including one additional constant to explain the size effect of micro-structures is introduced by Aifantis (1992). Although, in his study, stress concentration was quite removed at the tip of the crack and at the lines of dislocations/disclinations, the treatment of high-order derivatives of gradients in the governing equations was still needed. This is because the higher order derivatives are required ${C}^{1}$ continuous approximation space within the finite element framework. Later, Ru and Aifantis (1993) introduced the theorem to avoid continuity. In their approach, the fourth-order differential equations are parted into two parts of second-order differential equations. A staggered approach successfully handled boundary value problems of gradient elasticity in the framework of the finite element method (FEM) (Askes *et al*., 2008).

In this study, instead of FEM, we employ the smoothed finite element method (S-FEM) in the staggered approach for gradient elasticity. The smoothed finite element framework were introduced to improve the quality of the finite element solution obtained when using simplex elements, in particular, on triangles, the S-FEM variants such as the node-based S-FEM (NS-FEM) and edge-based S-FEM (ES-FEM) in 2D and NS-FEM and face-based S-FEM (FS-FEM) in 3D have shown to yield more accurate results than the traditional FEM. While cell-based S-FEM (CS-FEM) theoretically is similar to the traditional FEM and so does not relatively improve the solution over the triangular elements. These merits do come at a cost. The bandwidth of the the bilinear form computed using the ES-FEM/NS-FEM/FS-FEM is larger than the traditional FEM, which in turn increases the storage. However, there is a tradeoff between the desired accuracy and computational cost.

Since Liu *et al*. (2007a) introduced S-FEM, it has shown its toleration to shear/volumetric locking and heavily distorted meshes (Francis *et al*., 2022; Jiang *et al*., 2015; Lee *et al*., 2017). In addition, S-FEM has solved various engineering problems, in particular, discontinuity (Kshirsagar *et al*., 2021; Lee *et al*., 2023a; Surendran *et al*., 2021). Besides such aspects, S-FEM in general shows fast convergence and yields high accuracy than those of FEM in classical linear elasticity (Bordas and Natarajan, 2010; Liu *et al*., 2007b). Hence, S-FEM is suitable as a possible candidate for gradient elasticity to draw the size effect for micro-structures. Amongst different strain smoothing schemes, cell-based and edge-based strain smoothing methods are considered in this study.

This paper is organized as follows: the fundamental background of the stress-based gradient elasticity in the framework of finite element approximation is briefly introduced in Section 2 and the overview of smoothed finite element method is revisited in Section 3. Systematic numerical studies are conducted to draw the effect of the internal length scales on gradient elasticity in the context of the proposed S-FEM. Lastly, overall remarkable conclusions are discussed in the last section.

^{} 2. Theoretical Background of Gradient Elasticity

In this section, the governing equation and the discretized weak form of stress-based gradient elasticity are briefly discussed. Aifantis and his co-workers (Aifantis, 2003; Altan and Aifantis, 1997) suggested the following simple strong form of gradient elasticity:

where 𝜎 is the Cauchy stress, the infinitesimal strain is represented by $\u03f5=\frac{1}{2}(\nabla u+\nabla {u}^{T})$, $\mathrm{\u2102}$ is the elastic modulus, ${l}_{0}$ is a material length parameter/internal length scale which relates to micro-structural characteristics and $u$ is the displacement field. The equilibrium equations for a homogeneous body in the absence of inertia and body force is given by: $\nabla \sigma =0$. Then the fourth-order system of equation is derived as:

##### (2)

$\frac{1}{2}\mathrm{\u2102}[{\nabla}^{2}u+{\nabla}^{2}{u}^{T}-{l}_{0}^{2}{\nabla}^{2}({\nabla}^{2}u+{\nabla}^{2}{u}^{T}\left)\right]$Note that to solve Eq. (2) numerical technique needs ${C}^{1}$ continuous shape functions. However, by virtue of Ru-Aifiantis theorem, Eq. (2) can be rewritten in such ${C}^{0}$ continuous shape functions are required. The corresponding weak form is:

##### (3)

${\int}_{\Omega}\u03f5\left(v{)}^{T}\sigma \right(u)d\Omega ={\int}_{\Gamma}{v}^{T}\hat{t}d\Gamma $where $u$ is the trial function, $v$ is the test function and $\hat{t}$ is the traction. By substituting Eq. (1) to Eq. (3), we can obtain:

##### (4)

${\int}_{\Omega}\u03f5\left(v{)}^{T}\mathrm{\u2102}\right(1-{l}_{0}^{2}{\Omega}^{2}\left)\u03f5\right(u)d\Omega ={\int}_{\Gamma}{v}^{T}\hat{t}d\Gamma $Then, we can derive Eq. (4) by applying integration by part on the higher-order term to get the following equation:

##### (5)

${\int}_{\Omega}\u03f5\left(v{)}^{T}\mathrm{\u2102}\u03f5\right(u)d\Omega +{\int}_{\Omega}\nabla \u03f5(v{)}^{T}{l}_{0}^{2}\mathrm{\u2102}\nabla \u03f5(u)d\Omega =\phantom{\rule{0ex}{0ex}}{\int}_{\Gamma}{v}^{T}\hat{t}d\Gamma +{\int}_{\Gamma}\u03f5(v{)}^{T}{l}_{0}^{2}\mathrm{\u2102}\u03f5\left(u\right)d\Omega $Note that in a ${C}^{1}$ continuous formulation, it is typical to assume that this higher-order boundary condition vanishes at the boundary (Askes and Aifantis, 2002). However, in this study, due to the use of operator split, the higher-order boundary condition also vanishes as a result of the chosen boundary conditions. Then, Eq. (5) leads to $Ku=f$ (the system of linear equations), where:

##### (6a)

$K=\sum _{e}K{\phantom{\rule{.5em}{0ex}}}^{e}=\sum _{e}\left[{\int}_{{\Omega}^{h}}{B}^{T}\mathrm{\u2102}Bd\Omega +{\int}_{{\Omega}^{h}}\nabla {B}^{T}{l}_{0}^{2}\mathrm{\u2102}\nabla Bd\Omega \right]$where $B$ is the strain-displacement matrix. In the current approach, the strain-displacement matrix $B$ is built over smoothing domains, while traditionally it is computed on elements. The detailed computation of the $B$ matrix will be discussed in Section 3.

To solve Eq. (6), Ru and Aifantis (1993) introduced the operator split approach. In the context of finite element implementation, this approach includes two steps: computation of local displacements and computation of non-local stress fields. Firstly, the computation of local displacement fields is the same as classical elasticity as follows:

Once, the local displacement field is obtained, the non-local stress is calculated as:

where the non-local stress ${\sigma}^{g}$ is a nodal unknown field and $L$ is the linear differential operator given as in 2D:

##### (9)

$L=\left[\begin{array}{cc}\frac{\mathrm{\partial}}{\mathrm{\partial}x}& 0\\ 0& \frac{\mathrm{\partial}}{\mathrm{\partial}y}\\ \frac{\mathrm{\partial}}{\mathrm{\partial}y}& \frac{\mathrm{\partial}}{\mathrm{\partial}x}\end{array}\right]$The weak form of Eq. (8) can be given by:

##### (10)

${\int}_{\Omega}({v}^{T}{\sigma}^{g}+\nabla {v}^{T}{l}_{0}^{2}\nabla {\sigma}^{g})d\Omega ={\int}_{\Omega}{v}^{T}\mathrm{\u2102}Lud\Omega $where ${\sigma}^{g}=\sum {N}_{I}{s}_{I}^{g}$ with a set of shape functions $N$ as given in Eq. (11) and the nodal non-local stress ${s}_{I}^{g}$.

##### (11)

$N=\left[\begin{array}{ccc}{N}_{I}& 0& 0\\ 0& {N}_{I}& 0\\ 0& 0& {N}_{I}\end{array}\right],\phantom{\rule{.5em}{0ex}}I=1,\phantom{\rule{.5em}{0ex}}\dots ,\phantom{\rule{.5em}{0ex}}n$where $n$ is the number of nodes of the element. Finally, Eq. (10) can be rewritten using Eq. (1) and the non-local stress definition as follows:

##### (12)

$\left[{\mathit{\int}}_{\mathit{\Omega}}\mathit{(}{\mathit{N}}^{\mathit{T}}\mathit{N}\mathit{+}{\mathit{l}}_{\mathit{0}}^{\mathit{2}}\frac{\mathit{\partial}{\mathit{N}}^{\mathit{T}}}{\mathit{\partial}\hat{\mathit{x}}}\frac{\mathit{\partial}\mathit{N}}{\mathit{\partial}\hat{\mathit{x}}}\mathit{)}d\Omega \right]\phantom{\rule{.5em}{0ex}}{s}^{g}={\int}_{\mathrm{\Omega}}{N}^{T}\mathrm{\u2102}Lud\mathrm{\Omega}$where a summation $\hat{x}=x,y$. In this paper, we utilise the smoothed finite element method, which will be discussed in Section 3, to solve the local displacement field and use the finite element method to get the non-local stress. Interested readers can refer to Bagni and Askes (2015) and Lee *et al*. (2023b) and references therein for further details.

^{} 3. Smoothed finite element method: overview

In this section, the basic idea behind the gradient smoothing approach in the framework of the finite element method is briefly revisited. Liu *et al*. (2007a) introduced the gradient smoothing approach in the framework of the finite element method, often called smoothed finite element method (S-FEM) to improve linear triangular and bilinear quadrilateral elements. The one of remarkable aspects of S-FEM is the sub-division of finite elements, so-called sub-cells. Then, sub-cells are reconstructed to smoothing domains by different schemes where gradients are smoothed over the smoothing domains. Hence, strains are continuous over the smoothing domains but are discontinuous across the boundaries of smoothing domains, not elements (Liu and Nguyen, 2016).

Fig. 1 illustrates two main schemes of S-FEM that how to construct the smoothing domains: cell-based and edge-based S-FEM. Fig. 1(a) explains the concept of cell-based S-FEM (CS-FEM) and Fig. 1(b) is about edge-based S-FEM (ES-FEM), respectively. As shown in Fig. 1, each 3-node triangular element is split into three triangular sub-cells. Then, in the cell-based scheme, for example, $\u2206234$ is divided into three sub-cells: $\u220623{O}_{2}$, $\u220634{O}_{2}$ and $\u220624{O}_{2}$. In CS-FEM, each sub-cell becomes the smoothing domain where strains will be smoothed. ES-FEM also needs to split triangular elements into three triangular sub-cells. However, the smoothing domain is built by the position of the target edge. For instance, as shown in Fig. 1(b), when the target edge $k$ that the sharing boundary of $\u2206123$ and $\u2206234$ is chosen, the smoothing domain is combined by a sub-cell $\u220623{O}_{1}$ of $\u2206123$ and a sub-cell $\u220623{O}_{2}$ of $\u2206234$. Note that, when the target edge is the boundary of the domain, only one sub-cell is used as the smoothing domain in ES-FEM since the edge does not have a neighboring finite element.

In S-FEM approximation, the following infinitesimal strain ${\stackrel{~}{\u03f5}}^{h}$ over the smoothing domain ${\Omega}_{k}$ is considered:

##### (13)

${\stackrel{~}{\u03f5}}^{h}\left({x}_{k}\right)={\int}_{{\Omega}_{k}}{\u03f5}^{h}\left(x\right)\Phi \left(x\right)d\Omega ,\phantom{\rule{.5em}{0ex}}\forall x\in {\Omega}_{k}$where a point ${x}_{k}$ is located on the boundaries of the smoothing domain and the weight function 𝛷 is given in Eqs. (14a) and (14b):

##### (14b)

$\mathrm{\Phi}\left(x\right)=\left\{\begin{array}{l}\begin{array}{ll}1/{A}_{k}& x\in {\mathrm{\Omega}}_{k}\end{array}\\ \begin{array}{ll}0& \phantom{\rule{.5em}{0ex}}\phantom{\rule{.5em}{0ex}}\phantom{\rule{.5em}{0ex}}\phantom{\rule{.5em}{0ex}}\phantom{\rule{.5em}{0ex}}\phantom{\rule{.5em}{0ex}}\phantom{\rule{.5em}{0ex}}x\notin {\mathrm{\Omega}}_{k}\end{array}\end{array}\right.$Hence, Eq. (13) can be rewritten by means of the divergence theorem as follows:

##### (15)

$\u03f5\left({x}_{k}\right)=\frac{1}{{A}_{k}}{\int}_{{\Omega}_{k}}\u03f5\left(x\right)\Phi \left(x\right)d\Omega =\frac{1}{{A}_{k}}{\int}_{{\Gamma}_{k}}n\left(x\right){u}^{h}\left(x\right)d\Omega $where ${A}_{k}$ is the area of the smoothing domain and $n\left(x\right)$ is the outward normal vector given at Gauss points located on the middle point of boundaries of the smoothing domain as shown in Fig. 1. In terms of the nodal displacements, Eq. (15) is rewritten as:

##### (16)

$\stackrel{~}{\u03f5}\left(x\right)=\sum _{I\in {G}_{k}}\stackrel{~}{{B}_{I}}\left(x\right){u}_{I}^{h}$where ${G}_{k}$ is a set of nodes. The 2D smoothed strain-displacement matrix from Eq. (16) is defined as follows:

##### (17)

${\stackrel{~}{B}}_{IR}=\left[\begin{array}{cc}{\stackrel{~}{B}}_{Ix}& 0\\ 0& {\stackrel{~}{B}}_{Iy}\\ {\stackrel{~}{B}}_{Iy}& {\stackrel{~}{B}}_{Ix}\end{array}\right]$where

##### (18)

${\stackrel{~}{B}}_{Ii}=\frac{1}{{A}_{k}}{\int}_{{\Gamma}_{k}}{N}_{I}\left(x\right){n}_{i}\left(x\right)d\Gamma $where ${N}_{I}$ is the shape functions and ${\Gamma}_{k}$ is the boundaries of the smoothing domain and Eq. (19) denotes the 2D matrix form of the outward normal vectors. Since the 2D matrix form of the outward normal vectors is equivalent to the differential operator (Eq. (9)), Eq. (17) can be utilized for the computation of the stiffness matrix in the same manner as in FEM. As shown in Fig. 1 and Eq. (18), the proposed S-FEM does not require an explicit form of shape functions. In other words, there is no isoparametric mapping. In addition, the interior Gauss integration is altered to line integration with one Gauss point at the middle of boundaries of the smoothing domain (please refer to Fig. 1). These salient features of S-FEM bring marked benefits that it can avoid over-estimation of the stiffness and is immune to highly distorted meshes.

The discrete system equations for S-FEM is given by:

##### (20)

$\begin{array}{rl}& \stackrel{~}{K}u=f\\ & \stackrel{~}{K}={\int}_{\mathrm{\Omega}}{\stackrel{~}{B}}^{T}\left(x\right)\stackrel{~}{\mathrm{\u2102}B}\left(x\right)d\mathrm{\Omega}=\sum _{k=1}^{n}{\int}_{{\mathrm{\Omega}}_{k}}{\stackrel{~}{B}}^{T}\left(x\right)\mathrm{\u2102}B\left(x\right)d\mathrm{\Omega}\\ & =\sum _{k=1}^{n}{\stackrel{~}{B}}^{T}\left(x\right)\stackrel{~}{\mathrm{\u2102}B}\left(x\right){A}_{k}\end{array}$where $n$ is the number of target cells/edges for each scheme. The external force vector $f$ is obtained in the same manner in finite element approximation. In the present approximation, Eq. (7) is replaced by Eq. (20) for the local displacement fields. Eq. (12) is applied to each finite element in the current framework to obtain the non-local stress field.

^{} 4. Numerical tests

This section presents a series of two-dimensional tests with the proposed CS-FEM and ES-FEM. Fig. 2 illustrates the main process of the present framework. After generating triangular finite meshes, the smoothed stiffness matrix is constructed on either cell-based smoothing domains or edge-based smoothing domains. The non-local stress field is then calculated element-wise using Eq. (12).

A cantilever beam subject to end shear is validated to show the accuracy of the proposed S-FEM, compared to the analytical solution. Then, to capture the size effect at the dislocation lines and crack tip, a plate with a U-notched hole under far-field tension and a plate with a crack are examined. For all examples, 3-node triangular finite element meshes are used. To discuss the results, we use the following conventions:

∙ CS-FEM: cell-based smoothed finite element method

∙ ES-FEM: edge-based smoothed finite element method

∙ T3: 3-node triangular element

∙ DOFs: degrees of freedom

To examine tests, the proposed framework has been coded in MATLAB 2023b and the implementation was conducted on an Apple M2 chip with a clock speed of 3.49 GHz and 16 GB of memory.

4.1 Cantilever beam

This section tests a cantilever beam as the validation for the proposed smoothed FEM. Fig. 3 depicts the geometry of the beam and boundary conditions. The length of the beam is $L\in 0,5$ m and the height is given as $D\in -1,1$ m. The boundary conditions are imposed on the left end of the beam and the load $P$=1 N carries on the right side of the beam. Young’s modulus and Poisson’s ratio are $E$=1000 N/m^{2} and 𝜐=0.3, respectively. Note that the internal length scale is ${l}_{0}$=0.0. The exact displacements and stresses can be found in Augarde and Deeks (2008): ${\sigma}_{xx}=\frac{P(L-x)y}{I},\phantom{\rule{.5em}{0ex}}{\sigma}_{yy}=0,\phantom{\rule{.5em}{0ex}}{\tau}_{xy}=-\frac{P}{2I}\left(\frac{{D}^{2}}{4}-{y}^{2}\right),\phantom{\rule{.5em}{0ex}}{u}_{x}=-\frac{Py}{6EI}\left[(6L-3x)x+(2+\upsilon )\left({y}^{2}-\frac{{D}^{2}}{4}\right)\right]$ and ${u}_{y}=-\frac{P}{6EI}\left[3\upsilon {y}^{2}\right(L-x)+(4+5\upsilon )\frac{{D}^{2}x}{4}+(3L-x\left){x}^{2}\right]$. For this test, the following 9 levels of refined T3 elements are used: 306, 650, 1122, 1722, 2450, 3306, 4290, 5402 and 6642 DOFs. Fig. 4 shows the convergence of displacements and stresses for FEM, CS-FEM and ES-FEM. As expected, CS-FEM provides marginally better results in displacements and stresses to FEM; while ES-FEM yields more accurate results than other approaches.

4.2 Plate with a U-notched hole

Next, the infinite plate with a U-notched hole with the lateral extension is considered. As shown in Fig. 5, the length and height of the plate are $L$=5 m with the height of notch $h$=0.1 m and the radius of the hole is $r$=0.1 m, respectively. Lateral tension ${\sigma}_{xx}={10}^{3}$ N/m^{2} is implemented on the right edge of the plate. The same material properties used in Section 4.1 are also used for this test.

It is known that the infinite plate with a hole shows the stress concentration at the top of the hole because of a geometric discontinuity. Hence to quantify the effect of the stress concentration by the internal length scale, we measure the stress concentration factor ${K}_{t}$ given as (Pilkey and Pilkey, 2008):

${K}_{t}=\frac{{\sigma}_{max}}{{\sigma}_{nom}}={C}_{1+}{C}_{2}\left(\frac{h}{t}\right)+{C}_{3}\left(\frac{h}{t}{)}^{2+}{C}_{4}\right(\frac{h}{t}{)}^{3}$

where for $0.1\le h/r\le 2.0$, ${C}_{1}=0.955+2.169\sqrt{h/r}-0.081h/r,\phantom{\rule{.5em}{0ex}}{C}_{2}=-1.557-4.046\sqrt{h/r}+1.032h/r,\phantom{\rule{.5em}{0ex}}{C}_{3}=4.013+0.424\sqrt{h/r}-0.748h/r$ and ${C}_{4}=-2.461+1.538\sqrt{h/r}-0.236h/r$. ${\sigma}_{max}$ is obtained at the top of the notch where $y=h+r$. Therefore, the stress concentration factor of the plate is ${K}_{t}$=3.656.

Fig. 6 shows the convergence of the stress concentration factors with respect to the length of the plate and the internal length scales, ${l}_{0}$={0.001,0.005,0.01,0.05,0.1,0.5,1.0}. As shown in Fig. 6, ES-FEM gives a more accurate ${K}_{t}$ than CS-FEM. Furthermore, ${K}_{t}$ of both methods is decreased when the internal length scales are close to 1.0.

Fig. 7 illustrates stress ${\sigma}_{xx}$ distribution with the different length scales, ${l}_{0}$=0.0 and ${l}_{0}$=0.5. Figs. 7(a) and 7(b) show the stress distribution when ${l}_{0}$=0.0, i.e. the classical elasticity. Figs. 7(c) and 7(d) show the distribution of stress for CS-FEM and ES-FEM when ${l}_{0}$=0.5. It can be found that when the internal length scale is increased, the stress concentration is vanished and the values of stress are dramatically decreased.

4.3 Plate with a crack

Lastly, to verify the relaxation of stress concentration in stress-based gradient elasticity, we consider a plate with a crack, as shown in Fig. 8. The plate is a rectangular shape with dimensions of $L$=1 mm and the bottom edge is fully constrained. The crack implies as $d$=0.3 mm at $y$=1.0 mm. Two different material properties are assigned to two regions: (a) the left half ($x$=0~0.5 mm) has Young’s modulus ${E}_{1}$=17 GPa and Poisson’s ratio ${v}_{1}$=0.15 and (b) the right half ($x$=0.5~1.0 mm) has Young’s modulus ${E}_{2}$=200 GPa and Poisson’s ratio ${v}_{2}$=0.3, respectively. For this study, 8 different levels of the internal length scales are evaluated: ${l}_{0}$={0.0,0.001,0.005,0.01,0.05,0.1,0.5,1.0}. Both CS-FEM and ES-FEM simulations use 25600 elements (26130 DOFs).

Fig. 9 presents the stress distributions obtained using CS-FEM and ES-FEM for different internal length scales (51,842 degrees of freedom). As expected, the stress is highly concentrated at the crack tip and the bi-material interface. However, with increasing internal length scale, the stress becomes more evenly distributed throughout the domain, demonstrating the smoothing effect of the internal length scale in the proposed gradient elasticity.

Fig. 10 depicts the stress concentration reduction along the midlines of the plate ($y$=1.0 mm) for CS-FEM and ES-FEM. As expected, high stress concentrations are initially observed at the crack tip and the material interface, as shown in Fig. 10. However, the effect of the internal length scale becomes evident as the stress concentration alleviates. Notably, when the internal length scale approaches ${l}_{0}$=0.5, the stress distribution becomes finite, indicating a complete absence of sharp stress gradient.

^{} 5. Conclusions

This paper introduces CS-FEM and ES-FEM for two-dimensional stress-based gradient elasticity. The effectiveness of the proposed framework is demonstrated through various standard problems, consistently yielding accurate results. The influence of the internal length scale on stress concentration is investigated. The versatility of the framework is showcased by its applicability to problems with both smoothed and singular solutions. In particular, a single crack problem with bi-material properties is analyzed. Here, the internal length scale effectively mitigates the singularity typically observed at the crack tip and material interface. The proposed methods offer several advantages:

∙ avoids the need for explicit shape functions, thus, simplifying implementation;

∙ staggered scheme suits local displacement and non-local stress field computations, enhancing computational efficiency;

∙ implementation requires any further human intervention into the existing FEM code for gradient elasticity framework.

However, the proposed approach generally yields a higher computational cost than the standard FEM when employing the same number of DOFs. This increased cost derives from the subdivision of each finite element into sub-domains, which are utilized to construct smoothing domains. Consequently, the proposed S-FEM exhibits a wider stiffness bandwidth than FEM.

In future work, the application of gradient elasticity to functionally graded materials within the context of stress-constrained structural topology optimization will be explored.