## 2 Theory and Formulation
#### Nonlinear Temperature Changes
\begin{array}{ll}C_{ij}^{(r)}=\tau\Big{[}C_{ij}^{(r-1)}\times C_{ij}^{(1)}-C_{ ij}^{(r-1)}\big{/}\big{(}x_{i}-x_{j}\big{)}\Big{]}&i,j=1:n,\;\;i\
eq j\;\;\;\text{and}\;\;\;2\leq r \leq n-1\\ C_{ii}^{(r)}=-\sum\limits_{j=1,i\
eq j}^{n}C_{ij}^{(r)}&i,j=1:n\;\;\;\text{and} \;\;\;1\leq r\leq n-1\end{array}  \tag{18}\]  
In the current work, through x and y directions, a nonuniform set of points is presented as below [55]:  
\[ \begin{array}{ll}x_{i}=0.5\times a\times\Big{(}1-\cos\Big{(}\frac{(i-1)\pi} {(N_{i}-1)}\Big{)}\Big{)}&i=1:N_{i},\\ y_{i}=0.5\times b\times\Big{(}1-\cos\Big{(}\frac{(i-1)\pi}{(N_{i}-1)}\Big{)} \Big{)}&j=1:N_{j}\end{array}  \tag{19}\]  
Freedom degrees are represented as:  
\[u=U\times e^{i\omega t},v=V\times e^{i\omega t},w=W\times e^{i\omega t},\psi_ {x}=\Psi_{x}\times e^{i\omega t},\psi_{\theta}=\Psi_{\theta}\times e^{i\omega t}. \tag{20}\]  
Substitution of Equations (15) and (20) into governing equations, results in the following for natural frequency:  
\[\left\{\left[ \begin{array}{cc}[K_{dd}]&[K_{db}]\\ [K_{bd}]&[K_{bb}]\end{array} \right]\right.-\left[ \begin{array}{cc}[M_{dd}]&[ M_{db}]\\ [M_{bd}]&[M_{bb}]\end{array} \right]\omega^{2}\right\}\left\{ \begin{array}{c} \delta_{d}\\ \delta_{b}\end{array} \right\}=0 \tag{21}\]  
where, boundary and domain points are defined by subscripts \(b\) and \(d\), respectively. Additionally, the displacement vector is determined by the parameter of \(\delta\). It should be noted that in the stiffness matrix, all terms of the thermal loading in the previous section are added. The standard eigenvalue problem form of relation Equation (22) is stated as follows:  
\[ \begin{array}{l}[K^{*}][\delta_{i}]=\big{(}\omega^{2}\big{)}[M^{*}]\{\delta_ {i}\}\\ [K^{*}]=[K_{dd}-K_{db}K_{bb}{}^{-1}K_{bd}]\\ [M^{*}]=[M_{dd}-M_{db}K_{bb}{}^{-1}K_{bd}]\end{array}  \tag{22}\]  
It should be noted that by considering \(M^{*}=0\), we can obtain the critical temperature of the structure. It is noted that, in the critical temperature, thermal buckling occurs. For higher temperatures than the critical temperature, the structure encounters instability conditions, and that for the modeling of the structure, all designers should pay close attention to lower temperatures than the critical one. **Convergence Study**  
An adequate number of grid points is essential to attain correct results in the GDQ method in order to achieve convergence results. Different effective parameters of the system, such as boundary conditions (B.Cs) and materials, have been studied and compared. Furthermore, it could be observed that the system with four clamped edges has more stiffness than the system with (C-F-F-F) boundary conditions. Also, by adding the GPL reinforcing nanofillers to the system, the dimensionless natural frequency enhances. This means the GPLRC cylindrical panel becomes stiffer. Based on Table 1, the results are converged in 13 grid points. As Table 2 shows, the critical temperature results are converged in the 15 grid points. Furthermore, it should be expressed that the changes in the nondimensional length scale parameter have a direct influence on the critical temperature load. It means an increment in the nondimensional length scale leads to an increment trend in the critical temperature load. It should be noted that the dimensionless form of frequency and critical temperature (\(\Delta T^{*}_{Cr}\)) are presented as follows:  
\[\overline{\omega}=\omega h\sqrt{\rho/E}\;,\;\;\;\;\;\Delta T^{*}_{Cr}=100\times \Delta T_{cr}\times\overline{\alpha} \tag{23}\]  
For further validation and shape mode analysis of the GPLRC panel, finite element analyses were presented with the aid of ABAQUS, where shell elements (S4R) with a 4-node doubly curved panel, reduced integration, and finite membrane strains were used to create the mesh for the panel model. In addition, the perfect bonding between neighboring layers was considered. Furthermore, boundary conditions were applied to the nodes at the four edges of the GPLRC panel.