## 4 Numerical Experiments  
The enriched NEM module was developed in Fortran and combined into the previously developed PG-NEM code [30] having the stress recovery module [36]. To demonstrate and validate the present method, the crack problem of an infinite plate with collinear cracks depicted in Figure 5a is considered. Differing from the cracks in plates of finite size which are of great practical interest, the present infinite plate with periodic cracks provides the closed form solution. Meanwhile, for the numerical study, one may take two kinds of periodic finite crack problems. One is obtained when the plate is cut along EF and GH, while the other is obtained if the plate is cut along AB and CD. Since the first one has been widely taken for the numerical studies (i.e., [PERSON] and [PERSON] [37] and [PERSON] [38]), the second periodic analysis model is taken for the current numerical experiments. Figure 5b represents the detailed second periodic model, as well as a NEM grid having higher grid density toward the crack tip for an upper left darkened quarter. The SIF for this infinite plate with collinear cracks is expressed in the following closed form (the correction factor C of 1.02) [39]:  
\[K_{I}=\sigma_{\infty}\sqrt{\pi a}\left(\frac{2b}{\pi a}tan\frac{\pi a}{2b} \right)\cdot C \tag{43}\]  
with \(a\) being the half crack length, where a circular path \(\Gamma\) around the left crack tip has the relative radius \(r/a=0.0139\), and the near-tip stress distributions are to be evaluated along this circular path. The Young's modulus \(E\) is 200 GPa and Poisson's ratio \(\
u\) is 0.3. The problem was solved with 13 Gauss integration points, and the stress field approximated by PG-NEM was smoothened by the stress recovery technique [36] with the Laplace interpolation functions. The stress distributions were evaluated along the circular path \(\Gamma\) and are compared in Figure 6. It can be observed that the case without enrichment shows similar stress distributions to the exact ones, but it provides a significantly low stress level. Thus, one can realize that only the fine gradient grid cannot sufficiently capture the near-tip stress distribution. Meanwhile, it is clear that the case with enrichment remarkably improves the stress level. Therefore, it has been justified that the enrichment of interpolation functions is effective in enhancing the near-tip stress approximation. Figure 5: (**a**) Infinite plate with collinear cracks in plane strain state [2], (**b**) periodic model with double edge cracks and a gradient natural element method (NEM) grid (unit: m, \(N=2516\)). Figure 7a shows a rectangular FGM plate with an edge crack which is in the plane strain condition, where \(W\) and \(H/W\) are 1.0 and 8.0 units, while \(a/W\) is set variable for the parametric investigation as follows: 0.2, 0.3, 0.4, 0.5, and 0.6. The Young's modulus \(E\) varies along the \(x-\)axis in a form of  
\[E(x)=E_{1}exp(\beta x),\quad 0\leq x\leq W \tag{44}\]  
but [PERSON]'s ratio is uniform as 0.3. When \(E_{1}\) and \(E_{2}\) are specified to the elastic moduli at \(x=0\) and \(x=L\), respectively, the parameter \(\beta\) is calculated according to \(\beta=ln(E_{2}/E_{1})\), with \(E_{1}\) being 1.0 unit. Two material variation cases of \(E_{2}/E_{1}\) are taken as follows: \(E_{2}/E_{1}\) = 0.1 and 10.0, and two loading cases of tension in Figure 7b and bending in Figure 7c are considered. This example was initially studied by [PERSON] and [PERSON] [PERSON] and [PERSON] (2009), who also provided an analytical solution. The external loading is set by as follows: for tension and \(\sigma_{yy}(x,\pm 4)=\pm 1.0\) for tension and \(\sigma_{yy}(x,\pm 4)=\pm(-2x+1)\) for bending. Figure 6: Comparison of the stress distributions along the path \(\Gamma\) (\(r/a=0.0139\)): (**a**) \(\sigma_{xx}\), (**b**) \(\sigma_{yy}\), and (**c**) \(\tau_{xy}\). An upper grayed half is taken for the crack analysis from the problem symmetry, and the following symmetric constraint is specified to the symmetric line.