## 4 Natural Element Approximation of Thermal Buckling and Postbuckling  
The thermoelastic large deformation of an FG-GPLRC plate with a central crack is governed by the virtual work principle given by  
\[\delta U=0 \tag{27}\]  
In the PFF, the virtual strain energy \(\delta U\) of the cracked FG-GPLRC plate is calculated by [47]  
\[\delta U=\int_{-h/2+e}^{h/2-e}\int_{\mathcal{O}}{(1-\psi)^{2}\delta e^{T} \sigma\,dA\,dz} \tag{28}\]  
with \(e\) being the vertical distance to the neutral surface. Then, substituting Equations (22)-(26) into Equation (27) leads to the nonlinear simultaneous equations given by  
\[\left[\left(\mathbf{K}_{L,\sigma}+\sum_{e=1}^{M}\mathbf{K}_{L,s}^{e}\right)+\mathbf{K}_{ NL}\right]\overline{\mathbf{d}}=\mathbf{F} \tag{29}\]  
to compute the extended nodal vector \(\overline{\mathbf{d}}=\left\{\overline{b}_{1},\overline{b}_{2},\cdot\cdot\cdot, \overline{b}_{N},\overline{s}_{1}^{1},\cdot\cdot,\overline{s}_{2}^{4}\right\} ^{T}\) which is composed of the non-singular and singular nodal values. The linear and nonlinear stiffness matrices \(\mathbf{K}_{L,\sigma}\), \(\mathbf{K}_{L,s}^{e},\mathbf{K}_{NL}\) and the load vector \(\mathbf{F}\) will be defined later. Meanwhile, the thermal buckling and postbuckling of the FG-GPLRC plate under thermal loading is governed by the variation in total potential energy:  
\[\delta U-\delta W_{g}=0 \tag{30}\]  
where the variation in total potential energy \(\delta W_{g}\) due to the in-plane thermal loading is defined by  
\[\delta W=\int_{-h/2+e}^{h/2-e}\int_{\mathcal{O}}{(1-\psi)^{2}\delta\{u_{a, \beta}\}^{T}\mathbf{S}_{b}\,\delta\{u_{a,\beta}\}\,dA\,dz} \tag{31}\]  
with the \((9\times 9)\) matrix containing the thermal stresses  
\[\mathbf{S}_{b}= \begin{bmatrix}\mathbf{S}_{0}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{S}_{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&\mathbf{S}_{0}\end{bmatrix} \!,\,\mathbf{S}_{0}= \begin{bmatrix}\sigma_{ xx}&\sigma_{xy}&\tau_{xz}\\ \sigma_{yx}&\sigma_{yy}&\tau_{yx}\\ \tau_{xz}&\tau_{zy}&0\end{bmatrix}  \tag{32}\]  
and the displacement gradient vector \(\{u_{a,\beta}\}\) given by  
\[ \begin{array}{ll}\{u_{a,\beta}\}&=\begin{bmatrix}u_{x}\,u_{y}\,u_{z}\,v_{x} \,v_{gy}\,v_{z}\,w_{x}\,w_{y}\,w_{z}\end{bmatrix}^{T}\\ &=\mathbf{G}_{\mathbf{b}}\mathbf{b}_{I}+\hat{\mathbf{G}}_{I}^{k}\mathbf{s}_{I}^{k}\end{array}  \tag{33}\]In which \(\mathbf{G_{I}}\) and \(\mathbf{G_{I}^{k}}\) are the \((9\times 5)\) matrix containing the partial derivatives of L/I functions \(\varphi_{J}\) and \(\Phi_{I}^{k}\) respectively, and which are defined by  
\[\mathbf{G_{I}}= \begin{bmatrix}\hat{G_{I1}}&0&0&\mathbf{G_{I2}}&0\\ 0&\mathbf{G_{I1}}&0&0&\mathbf{G_{I2}}\\ 0&0&\mathbf{G_{I1}}&0&0\end{bmatrix} ,\ \hat{\mathbf{G}_{I}}^{k}= \begin{bmatrix}\hat{\mathbf{G}_{I1} }^{k}&0&0&\hat{\mathbf{G}_{I2}}^{k}&0\\ 0&\hat{\mathbf{G}_{I1}}^{k}&0&0&\hat{\mathbf{G}_{I2}}^{k}\\ 0&0&\hat{\mathbf{G}_{I1}}^{k}&0&0\end{bmatrix}  \tag{34}\]  
with  
\[\mathbf{G_{I1}}=\left\{\frac{\partial\varphi_{I}}{\partial x}\,\frac{\partial \varphi_{I}}{\partial y}\,0\right\}^{T},\mathbf{G_{I2}}=\left\{\frac{2\,\partial \varphi_{I}}{\partial x}\,z\,\frac{\partial\varphi_{I}}{\partial y}\,\varphi_ {I}\right\}^{T} \tag{35}\]  
Here, \(\hat{\mathbf{G}_{I}^{k}}\) and \(\hat{\mathbf{G}_{I2}^{k}}\) are defined by replacing \(\varphi_{I}\) in Equation (35) with \(\Phi_{I}^{k}\). The thermal stresses in Equation (32) are calculated using Equations (25) and (26) after the thermoelastic analysis in Equation (29). Meanwhile, in this study, the phase field \(\psi(\mathbf{x})\) is also approximated using the same L/I functions taken for the displacement approximation. However, the crack line may not coincide with any one in vertical, horizontal or inclined lines within the NEM grid, as illustrated in Figure 4, where a red line indicates an inclined central crack. In this case, the phase field approximated using the NEM grid cannot represent such a crack line without any extra numerical technique.