## 5 Numerical Results
The iteration compensation factor \(\sigma\) was set to \(\sigma=\left(1+\overline{w}_{c}^{2}\right)\) with the calibrated central deflection \(\overline{w}_{c}=w_{c}/h\), and which was chosen by referring to the solutions of [PERSON] et al. [19] and [PERSON] et al. [54]. This factor was adopted because the nonlinear stiffness matrix \(\mathbf{K}_{NL}\) in Equation (41) is modified during the nonlinear iteration, which differs from the linear shear stiffness matrix \(\mathbf{K}_{L,s}^{\sigma}\) in Equation (40), which is unchanged during the iteration. The computed thermal postbuckling path is represented in Figure 6 with two reference solutions, where \(\Delta T_{cr,L}\) and \(\Delta T_{cr,NL}\) indicate the linear and nonlinear CBTR. [PERSON] et al. [54] obtained the path by combining a Galerkin differential quadrature iteration algorithm and [PERSON]'s higher-order SDT. \begin{table}
\begin{tabular}{c c c c c c c c} \hline \multirow{2}{*}{**Items**} & \multicolumn{6}{c}{**Grid Density (\(\mathbf{N\times N}\))**} \\ \cline{2-9}  & **11\(\times\)11** & **13\(\times\)13** & **15\(\times\)15** & **17\(\times\)17** & **19\(\times\)19** & **21\(\times\)21** & **23\(\times\)23** \\ \hline \(\lambda_{T}\) & 0.1705 & 0.1542 & 0.1435 & 0.1361 & 0.1307 & 0.1264 & 0.1229 \\ \(\lambda_{T}^{rel}\)(\(\%\)) & 38.731 & 25.478 & 16.762 & 10.740 & 6.347 & 2.848 & - \\ \hline \end{tabular}
\end{table} 
Table 1: Convergence of CBTP \(\lambda_{T}=10^{3}\alpha T_{cr}\) to the grid density for FIMM simply-supported intact isotropic plate under a uniform temperature rise (\(a/b=1\), \(a/h=100,\
u=0.3\)). \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Methods** & _a/h = 10_ & _a/h = 20_ & _a/h = 100_ \\ \hline Present & 11.92 & 3.106 & 0.1264 \\ [PERSON] et al. [19] & 11.98 & 3.119 & 0.1265 \\ [PERSON] et al. [52] & 11.99 & 3.123 & 0.1272 \\ [PERSON] et al. [53] & 11.83 & 3.089 & 0.1271 \\ \hline \end{tabular}
\end{table} 
Table 2: Comparison of CBTP \(\lambda_{T}\) for FIMM simply-supported intact isotropic plates subject to uniform temperature rise (\(a/b=1\), \(\
u=0.3\)). Table 3 compares the CBTRs of a simply-supported Al/ZrO\({}_{2}\) functionally graded square plate with an inclined central crack. The volume fraction \(V_{c}(z)\) of ZrO\({}_{2}\) is expressed as \(V_{c}(z)=\left(0.5+z/h\right)^{n}\), and the elastic properties are \(E_{m}=70\) GPa, \(v_{m}=0.3\), \(\alpha_{m}=23\times 10^{-6}/^{\circ}\)C for Al and \(E_{c}=151\) GPa, \(\
u_{c}=0.3\), \(\alpha_{c}=10\times 10^{-6}/^{\circ}\)C for ZrO\({}_{2}\), respectively. The sides and thickness of the plate are \(a=b=1.0\) m and \(a/h=10\), and the crack angle and relative length are \(\alpha=90^{\circ}\) and \(c/a=0.6\). The CBTR values were computed for five different ceramic indices, \(n=0\), 0.5, 1, 2, and 5, and compared with the reference solutions. The CBTR has been shown to decrease uniformly proportional to the ceramic volume fraction index \(n\) because the increase in the volume fraction of metals with a higher coefficient of thermal expansion decreases the thermal buckling stiffness of the plate as the index \(n\) increases [55]. Meanwhile, good agreement between the present method and two reference methods is clearly seen such that the maximum relative difference is 5.585% at \(n=0.5\). [PERSON] et al. [56] obtained solutions by applying the phase-field theory to first-order SDT, while [PERSON] et al. [57] combined the extended isogeometric analysis (XIGA) and level sets into first-order SDT. The verification results show that the present phase-field crack model can stably predict the thermal buckling behavior of composite plates containing cracks. Next, the CBTRs of intact and immovable simply-supported FG-GPLRC plates with a side length \(a=b=1.0\) m and the GPL mass fraction \(g_{GPL}^{*}=0.3\%\) are computed and compared in Table 4. The CBTRs were computed by changing the thickness ratio \(a/h\) and the GPL distribution pattern for PIMM and FIMM simply-supported BCs.