## 3 NE Approximation of Nonlinear Natural Vibration Using Phase-Field Crack Model  
Referring to Figure 4, the panel neutral surface \(\varpi\) is uniformly divided into a finite number of three-node Delaunay triangles. As mentioned above, the crack line is modeled using the phase field \(\phi(x)\) and the crack-induced singular field is represented by adding the crack-tip singular functions without refining the local NEM grid in the vicinity of crack tip. Thus, in an extended NEM, the displacement \(\mathbf{u}(x,y,z)\) is approximated as  
\[\mathbf{u}^{h}(\mathbf{x})=\sum_{J=1}^{N}\left(\mathbf{\eta}_{J}+z\,\mathbf{d}_{J}\right)\, \Psi_{J}(x,y)+\sum_{I=1}^{2}\sum_{k=1}^{4}\left(\hat{p}_{I}^{k}+z\,\hat{\mathbf{d}} _{I}^{k}\right) \begin{cases}Q_{I}^{k}(r,\theta)\\ Q_{I}^{k}(r,\theta)\\ R_{I}^{k}(r,\theta)\end{cases}  \tag{25}\]  
using the L/I functions \(\psi_{J}(x,y)\)[46; 47] and the crack-tip singular functions \(Q_{k}(r,\theta)\) and \(R_{k}(r,\theta)\)[38]. \[\left\{Q_{I}^{k}(r,\theta)\right\}=\left\{\sqrt{\tau}sin\left(\frac{\theta}{2 }\right),\,\sqrt{\tau}cos\left(\frac{\theta}{2}\right),\,\sqrt{\tau}sin\left( \frac{\theta}{2}\right)sin(\theta),\,\sqrt{\tau}cos\left(\frac{\theta}{2} \right)sin(\theta)\right\} \tag{26}\]  
\[\left\{R_{I}^{k}(r,\theta)\right\}=\left\{\sqrt{\tau}sin\left(\frac{\theta}{2 }\right),\,0,\,0,\,0\right\} \tag{27}\]  
The singular functions are defined using two distinct crack-tip polar coordinates \((r,\theta)_{I}(I=1,\,2)\), as depicted in Figure 3. Here, \(\mathbf{b}_{J}=\left(\mathbf{r}_{J},\,\mathbf{d}_{J}\right)\) are the non-singular nodal displacement vectors at node \(J\) within the NEM grid \(\Im_{C}\) generated with \(N\) nodes and Delaunay triangles. And, \(\mathbf{s}_{l}^{k}=\left(\mathbf{\hat{r}}_{l}^{k},\mathbf{\hat{d}}_{l}^{k}\right)\) indicates the singular vector corresponding to the k-th singular function. The derivation and manipulation of L/I functions on the curved cylindrical surfaces are troublesome, so the physical NEM grid \(\Im_{C}=[0,S]\times[0,L]\) on the panel neutral surface is transformed into 2-D rectangular NEM grid \(\Im_{R}=[0,\theta_{0}]\times[0,L]\) for the computation. Then, according to the geometry transformation \(T_{C}\) and the chain rule summarized in Appendix A, the enriched NE approximation of the bending-membrane strain \(\varepsilon\) in Equation (18) and the transverse shear (T/S) strain \(\gamma\) in Equation (19) ends up with  
\[\varepsilon^{h}=\sum_{J=1}^{N}\widetilde{\mathbf{L}}\phi_{J}\mathbf{b}_{J}+\sum_{l=1}^ {2}\sum_{k=1}^{4}L\,\Phi_{J}^{k}\mathbf{s}_{l}^{k}=\sum_{J=1}^{N}B_{J}\mathbf{b}_{J}+ \sum_{l=1}^{2}\sum_{k=1}^{4}H_{l}^{k}\mathbf{s}_{l}^{k} \tag{28}\]  
\[\gamma^{h}=\sum_{I=1}^{N}\widetilde{\mathbf{L}}_{s}\phi_{J}\mathbf{b}_{J}+\sum_{I=1}^ {2}\sum_{k=1}^{4}L_{s}\Phi_{I}^{k}\mathbf{s}_{I}^{k}=\sum_{J=1}^{N}B_{J}^{l}\mathbf{b} _{J}+\sum_{I=1}^{2}\sum_{k=1}^{4}H_{l}^{k}\mathbf{s}_{I}^{k} \tag{29}\]  
with \(\Phi_{I}^{k}=\left\{Q_{I}^{k},\,Q_{I}^{k},\,R_{I}^{k}\right\}\). Here, \(\mathbf{B}_{I}=\widetilde{\mathbf{L}}\phi_{J}\) is computed on the 2-D rectangular NEM grid while \(\mathbf{H}_{I}^{k}=\mathbf{L}\Phi_{I}^{k}\) is computed directly on the panel neutral surface (Similarly for \(\mathbf{B}_{s}^{J}\) and \(\mathbf{H}_{SI}^{k}\)). Note that the first term on the right-hand side in Equation (29), which is approximated with C\({}^{0}\)-L/I functions \(\phi_{J}\), may suffer from numerical locking when the deformation is bending-dominated [39; 40]. This problem can be effectively suppressed by indirectly interpolating this non-singular term of T/S strain according to the MITC3+ shell approach [41], as addressed in Appendix B.