## 2 Equation of Motion for Nuclear Rod Bundles Surrounded by Viscous Fluid
### Numerical Model of Nuclear Rod Bundles  
Flexural wave propagations were calculated using SEM to evaluate the structural instability of nuclear rod bundles in the fluid medium. The spectral element equation of the weak form based on a weighted-integral statement was derived. An eigenvalue problem was solved to determine the wavenumber of the flexural wave for the nuclear rod bundles in the fluid medium. The effects of added mass and fluid damping were investigated by changing the number of structures. The transfer function of the flexural vibration that was induced by incident pressure wave was calculated to evaluate the structural instability of the nuclear rod bundles. Figure 3: Definition diagram, coordinate systems and symbols used in the determination of the streamline and pressure of the surrounding fluid of the cylinders. Deriving the dynamic shape functions depended on the frequency domain, and the general solution of Equation (1) was assumed as  
\[w_{i}(x,t)=\beta_{i}Ce^{-\mathrm{i}kx}e^{\mathrm{i}\omega t}, \tag{12}\]  
where \(\beta\) is the proportional coefficient for each general solution and \(k\) is the wavenumber. Substituting Equation (1) into (12), the eigenvalue problem was obtained by  
\[(k^{4}E\mathbf{I}-\omega^{2}(\rho_{s}A_{s}\mathbf{I}+\mathbf{M}_{a}))\mathbf{ \beta}=0, \tag{13}\]  
where \(\mathbf{I}\) is the unit matrix of \(N\) by \(N\), \(\mathbf{\beta}=\left\{ \begin{array}{cccc}\beta_{1}&\beta_{2}&\cdots&\beta_{N }\end{array} \right\}^{T}\) and \(\mathbf{M}_{a}\) is the matrix of added masses as  
\[\mathbf{M}_{a}=\left[ \begin{array}{cccc}M_{1,1}&M_{1,2}&\cdots&M_{1,N}\\ M_{2,1}&M_{2,2}&&M_{2,N}\\ \vdots&&\ddots&\vdots\\ M_{N,1}&M_{N,2}&\cdots&M_{N,N}\end{array} \right].\]  
Equation (13) was rewritten in a linearized form by  
\[\mathbf{C}\mathbf{x}=k\mathbf{J}\mathbf{x}, \tag{14}\]  
where  
\[\mathbf{C}=\left[ \begin{array}{cccc}0&\mathbf{I}&0&0\\ 0&0&\mathbf{I}&0\\ 0&0&0&\mathbf{I}\\ \omega^{2}(\rho_{s}A_{s}\mathbf{I}+\mathbf{M}_{a})&0&0&0\end{array} \right],\]  
and  
\[\mathbf{J}=\left[ \begin{array}{cccc}\mathbf{I}&0&0&0\\ 0&\mathbf{I}&0&0\\ 0&0&\mathbf{I}&0\\ 0&0&0&E\mathbf{I}\end{array} \right].\]  
After the QR algorithm was applied to numerically solve Equation (14) using MATLAB, \(4N\) wavenumbers \(k_{i}\), proportional coefficients \(\beta_{i}\) determined by the material properties and frequency were obtained. While using the \(4N\) wavenumbers, the general solutions of Equation (1) were obtained as  
\[W_{i}(x)=\mathbf{e}(x,\omega)\mathbf{B}_{i}(\omega)\mathbf{a}, \tag{15}\]  
where  
\[\mathbf{a}=\left\{ \begin{array}{cccc}a_{1}&a_{2}&a_{3}&\cdots&a_{4N}\end{array}  \right\}^{T},\]  
and  
\[\mathbf{B}_{i}(\omega)=\left[ \begin{array}{cccc}\beta_{i,1}&&&&\\ &\beta_{i,2}&&\\ &&\ddots&\\ &&&\beta_{i,4N}\end{array} \right].\]The spectral nodal displacement of the multiple-beam element in the flexural direction is expressed by  
\[\mathbf{d}=\Big{\{} \begin{array}{cccccccc}W_{1}(0)&W^{\prime}{}_{1}(0)&\cdots&W_ {N}(0)&W^{\prime}{}_{N}(0)&W_{1}(L)&W^{\prime}{}_{1}(L)&\cdots&W_{N}(L)&W^{\prime} {}_{N}(L)\end{array} \Big{\}}^{T}. \tag{16}\]  
Substituting Equation (15) into (16),  
\[\mathbf{d}=\mathbf{H}(\omega)\mathbf{a}, \tag{17}\]  
where  
\[\mathbf{H}(\omega)=\left[