## 3 Natural Element Approximation of Strain and Stress Fields  
By adopting the first-order SDT, the displacement field \(\mathbf{u}=\left\{u,v,w\right\}^{T}\) of an FG-GPLRC plate subjected to thermal loading is expressed as  
\[\left\{ \begin{aligned} & u\\ & w\end{aligned} \right\}_{(x,y,z)}=\left\{ \begin{aligned} & u_{0}\\ & w_{0}\end{aligned} \right\}_{(x,y)}+z\cdot\left\{ \begin{aligned} & \beta_{x}\\ &\beta_{y}\\ & 0\end{aligned} \right\}_{(x,y)} \tag{11}\]  
Here, the nodal translation vector \(\mathbf{a}=\left\{u_{0},v_{0},w_{0}\right\}^{T}\) and the nodal rotation vector \(\mathbf{r}=\left\{\beta_{x},\beta_{y}\right\}^{T}\) are combined into a displacement nodal vector \(\mathbf{b}=\left(u_{0},v_{0},w_{0},\beta_{x},\beta_{y}\right)^{T}\) for description convenience. Figure 2: The phase field \(\psi(\mathbf{x})\) within the cracked FG-GPLRC plate. The thermally induced large deflection of the FG-GPLRC plate is represented by von Karman nonlinearity. Then, the compatibility relations between displacement and strain are expressed as  
\[\varepsilon=\left\{ \begin{array}{c}\frac{\partial u_{0}}{\partial x}+0.5\vec{w}_ {x,x}\frac{\partial w}{\partial x}\\ \frac{\partial w}{\partial y}+0.5\vec{w}_{y}\frac{\partial w}{\partial y}\\ \frac{\partial w_{0}}{\partial x}+\frac{\partial u_{0}}{\partial y}+0.5\Big{(} \vec{w}_{x,x}\frac{\partial w}{\partial x}+\vec{w}_{y}\frac{\partial w}{ \partial y}\Big{)}\\ \end{array} \right\}+z\cdot\left\{ \begin{array}{c}\frac{\partial\beta_{x}}{ \partial x}\\ \frac{\partial w_{y}}{\partial y}\\ \frac{\partial w_{x}}{\partial y}+\frac{\partial\beta_{y}}{\partial x}\\ \end{array} \right\}=(\vec{H}_{L}+\vec{H}_{NL})\vec{b} \tag{12}\]  
\[\gamma=\left\{ \begin{array}{c}\beta_{y}+\frac{\partial w_{0}}{\partial y} \\ \beta_{x}+\frac{\partial w_{0}}{\partial x}\end{array} \right\}=\vec{H}_{s}\vec{b} \tag{13}\]  
with \(\varepsilon=\left\{\varepsilon_{xx},\varepsilon_{yy},2\varepsilon_{xy} \right\}^{T}\) and \(\gamma=\left\{\gamma_{yz},\gamma_{zx}\right\}^{T}\). \(\vec{H}_{L}\), \(\vec{H}_{NL}\) and \(\vec{H}_{s}\) indicate the \((3\times 5)\) and \((2\times 5)\) partial derivative matrices defined by  
\[\vec{H}_{L}= \begin{bmatrix}H_{x}&0&0&z\cdot H_{x}&0\\ 0&H_{y}&0&0&z\cdot H_{y}\\ H_{y}&H_{x}&0&z\cdot H_{y}&z\cdot H_{x}\end{bmatrix}  \tag{14}\]  
\[\vec{H}_{NL}= \begin{bmatrix}0&0&\vec{w}_{x}H_{x}/2&0&0\\ 0&0&\vec{w}_{y}H_{y}/2&0&0\\ 0&0&\left(\vec{w}_{x}H_{y}+\vec{w}_{y}H_{x}\right)/2&0&0\end{bmatrix}  \tag{15}\]  
\[\vec{H}_{s}= \begin{bmatrix}0&0&H_{y}&0&1\\ 0&0&H_{x}&1&0\end{bmatrix}  \tag{16}\]  
with \(H_{x}=\partial/\partial x\) and \(H_{y}=\partial/\partial y\). Here, \(\vec{w}_{x}\) and \(\vec{w}_{y}\) are the partial derivatives of \(w\) which are given a priori, as will be explained later. Then, the strain-stress constitutive relations, considering the thermal effect, are expressed as  
\[\sigma=\frac{E}{1-\
u^{2}} \begin{bmatrix}1&\
u&0\\ \
u&1&0\\ 0&0&(1-\
u)/2\end{bmatrix}\begin{cases}\varepsilon_{xx}-\alpha_{t}\Delta T\\ \varepsilon_{yy}-\alpha_{t}\Delta T\\ 2\varepsilon_{xy}\end{cases} =\vec{D}[(\vec{H}_{L}+\vec{H}_{NL})\vec{b}- \alpha_{t}\Delta T] \tag{17}\]  
\[\tau= \begin{bmatrix}G&0\\ 0&G\end{bmatrix}\begin{cases}\gamma_{yz}\\ \gamma_{zx}\end{cases} =\vec{D}_{s}\vec{H}_{s}\vec{b} \tag{18}\]  
with \(\sigma=\left\{\sigma_{xx},\sigma_{yy},\sigma_{xy}\right\}^{T}\), \(\tau=\left\{\tau_{yz},\tau_{zx}\right\}^{T}\) and \(\vec{I}=\left\{1,1,0\right\}^{T}\). For the natural element approximation of the thermally induced displacement field, the neutral surface \(\omega\) is uniformly divided into three-node Delaunay triangles, as shown in Figure 3. Delaunay triangles and nodes generate an NEM grid, and Laplace interpolation (L/I) functions and their domains are represented, together with a central crack in red. The abrupt increase in the stress distribution, called the singular behavior, in the vicinity of the crack tips, is represented by adding the crack-tip singular functions instead of the local grid refinement near the crack tip.