## 3 NE Approximation of Nonlinear Natural Vibration Using Phase-Field Crack Model
The analytical calculation of Equations (A5) and (A6) in Appendix B using Equations (19) and (22), together with the chain rule between two coordinates \((x,y)\) and \((\xi,\eta)\) in Figure 4, results in  
\[\gamma_{\varepsilon}^{h}=\mathbf{\hat{B}}_{\varepsilon}\mathbf{b}_{\varepsilon} \tag{30}\]  
Here, \(\mathbf{\hat{B}}_{\varepsilon}\) is the \((2\times 15)\) triangle-wise matrices in function of \(\xi,\eta,z\), and \(R\), and \(\mathbf{b}_{\varepsilon}=\left\{\mathbf{b}_{1}^{\varepsilon},\mathbf{b}_{2}^{\varepsilon}, \mathbf{b}_{3}^{\varepsilon}\right\}\) are the \((15\times 1)\) non-singular triangle-wise nodal vectors. Next, the nonlinear vibration of the FG-GPLRC cylindrical panel is governed by the dynamic form of the energy principle given by  
\[\delta U+\delta\,T=0 \tag{31}\]  
Figure 4: A geometry transformation \(T_{C}\) between the curved and planar NEM grids and L/I functions \(\phi_{I}(\xi_{1},\xi_{2})\) defined on a 2-D rectangular plane. In accordance with the phase field concept, the virtual strain energy \(\delta U\) and the virtual kinetic energy \(\delta T\) of cracked cylindrical panel are defined by  
\[\delta\,U=\int_{-h/2+e}^{h/2-e}\int_{\alpha}{(1-\phi)^{2}\delta\varepsilon^{T} \sigma\,dAdz} \tag{32}\]  
\[\delta\,T=\int_{-h/2+e}^{h/2-e}\int_{\alpha}{(1-\phi)^{2}\delta d^{T}\mathbf{m} \mathbf{\bar{d}}\,dAdz} \tag{33}\]  
where \(e\) denotes the radial distance of the neutral surface from the panel mid-surface and \(\mathbf{m}\) is a \((5\times 5)\) symmetric matrix given by  
\[\mathbf{m}=\rho \begin{bmatrix}\mathbf{I}&\mathbf{m}_{1}^{T}\\ \mathbf{m}_{1}&\mathbf{m}_{2}\end{bmatrix} ,\;\mathbf{m}_{1}= \begin{bmatrix}z&0&0\\ 0&z&0\end{bmatrix}  \tag{34}\]  
with \(\mathbf{m}_{2}=diag\big{(}z^{2},z^{2}\big{)}\) and the identity matrix \(\mathbf{I}\). Figure 5 represents a NEM grid composed of Delaunay triangles and grid points, which will be described later in detail. In the present study, the phase field \(\phi(\mathbf{x})\) introduced in Figure 4 is approximated using L/I functions defined on the NEM grid. So, the crack line in red may not have lied on the vertical, inclined, or horizontal lines within the NEM grid, as depicted in Figure 5. In such a situation, an additional numerical technique is needed to effectively approximate the phase field. In this study, the crack coordinates \(\mathbf{x}\in\Re_{C}^{2}\) including the crack line transforms the virtual phase field coordinates \(\mathbf{\bar{x}}\in\Re_{P}^{2}\) such that the crack line lies on the inclined line within the NEM grid. \[\mathbf{x}\in\Re_{C}^{2}\;\rightarrow\;\mathbf{\bar{x}}\in\Re_{P}^{2} \tag{35}\]  
\[ \begin{cases}\mathbf{\bar{x}}\\ \mathbf{\bar{y}}\end{cases} = \begin{cases}x_{o}\\ y_{o}\end{cases} +\mathbf{\bar{R}}\cdot \begin{cases}cos(\mathbf{\bar{\theta}})\\ sin(\mathbf{\bar{\theta}})\end{cases} ,\;\mathbf{\bar{\theta}}=tan^{-1}\Big{(}\frac{y}{ x}\Big{)}-\big{(}\alpha+\theta_{G}\big{)} \tag{36}\]  
with \(\mathbf{\bar{R}}^{2}=(x-x_{0})^{2}+(y-y_{0})^{2}\). Next, substituting Equations (32) and (33) into Equation (31) through the constitutive relations (23) and (24), under the assumption harmonic motion \(\mathbf{d}=\mathbf{\bar{d}}\cdot e^{|\omega t}\), one obtains the nonlinear modal equation given by  
\[\left[\left(\mathbf{K}_{L,\sigma}+\sum_{e=1}^{M}\mathbf{K}_{L,s}^{\epsilon}\right)+ \mathbf{K}_{NL}\right]\mathbf{\bar{d}}\;-\omega^{2}\mathbf{M}\mathbf{\bar{d}}=\mathbf{0} \tag{37}\]  
Figure 5: Rotation of inclined central crack for the phase field approximation in NEM.