## 2 One-Dimensional Formulation
### Finite Element Implementation  
The finite element method is employed in this paper for the ease of computation. In consideration of the second-order partial derivative in Equation (11), the quadratic interpolation function is used to approximate the axial displacement field within an element.  
\[\rmx=\bf{N}\rm{\bf X}=\left[\begin{array}{cccc}\rm{\bf N}_{1}&\rm{\bf N}_{2}& \rm{\bf N}_{3}\end{array}\right]\rmX \tag{12}\]where \(\mathbf{N}\) and \(\mathbf{X}\) represent the shape function matrix and the nodal generalized displacement vector, respectively. The submatrices in Equation (12) are written as follows:  
\[\mathbf{N}_{1}=\bigg{(}\frac{2z^{2}}{l^{2}}-\frac{3z}{l}+1\bigg{)}\mathbf{I}_{N },\ \mathbf{N}_{2}=\bigg{(}\frac{4z}{l}-\frac{4z^{2}}{l^{2}}\bigg{)}\mathbf{I}_{N}, \ \mathbf{N}_{3}=\bigg{(}\frac{2z^{2}}{l^{2}}-\frac{z}{l}\bigg{)}\mathbf{I}_{N} \tag{13}\]  
where the coordinate \(z\) varies from 0 to \(l\) within an element, and \(\mathbf{I}_{N}\) denotes the \(N\times N\) identity matrix. The vector \(\mathbf{X}\) in Equation (12) is divided into three submatrices.  
\[\mathbf{X}=\bigg{[}\begin{array}{cccc}\mathbf{X}_{1}^{\mathbf{T}}&\mathbf{ X}_{2}^{\mathbf{T}}&\mathbf{X}_{3}^{\mathbf{T}}\end{array}\bigg{]}^{\mathsf{T}}, \ \mathbf{X}_{i}=\bigg{[}\begin{array}{cccc}X_{1(i)}&X_{2(i)}&\cdots&X_{N(i)} \end{array}\bigg{]}^{\mathsf{T}},\ i=1,\ 2,\ 3 \tag{14}\]  
where the subscripts 1, 2, and 3 indicate two end and one middle nodes of an element. Specially,  
\[\begin{array}{l}\int\limits_{A}\Big{[}\mathbf{x}^{\mathsf{T}}\mathbf{b}^{ \mathsf{T}}\mathbf{H}^{\mathsf{T}}\overline{\sigma}\Big{]}_{z=z_{1}}^{z=z_{2} }\mathbf{d}A=\int\limits_{A}\big{[}\mathbf{x}^{\mathsf{T}}\big{]}_{z=z_{2}} \mathbf{\psi}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}[\overline{\sigma}]_{z=z_{2} }\mathbf{d}A-\int\limits_{A}\big{[}\mathbf{x}^{\mathsf{T}}\big{]}_{z=z_{1}} \mathbf{\psi}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}[\overline{\sigma}]_{z=z_{1} }\mathbf{d}A\\ =\int\limits_{A}\mathbf{X}_{3}^{\mathsf{T}}\mathbf{\psi}^{\mathsf{T}}\mathbf{H }^{\mathsf{T}}[\overline{\sigma}]_{z=z_{2}}\mathbf{d}A-\int\limits_{A} \mathbf{X}_{1}^{\mathsf{T}}\mathbf{\psi}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}} [\overline{\sigma}]_{z=z_{1}}\mathbf{d}A=\int\limits_{A}\mathbf{X}^{\mathsf{ T}}\Big{[}\overline{\mathbf{N}}^{\mathsf{T}}\mathbf{\psi}^{\mathsf{T}}\mathbf{H}^{ \mathsf{T}}\overline{\sigma}\Big{]}_{z=z_{1}}^{z=z_{2}}\mathbf{d}A\end{array} \tag{15}\]  
where \(\Big{[}\overline{\mathbf{N}}^{\mathsf{T}}\Big{]}_{z=z_{2}}=[\mathbf{O}_{N}, \mathbf{O}_{N},\mathbf{I}_{N}]\), \(\Big{[}\overline{\mathbf{N}}^{\mathsf{T}}\Big{]}_{z=z_{1}}=[\mathbf{I}_{N}, \mathbf{O}_{N},\mathbf{O}_{N}]\) and \(\mathbf{O}_{N}\) denotes the \(N\times N\) null matrix.  
Substituting Equations (12) and (15) into Equation (11), and referring to the relationship below:  
\[\delta\mathbf{x}^{\mathsf{T}}\
eq 0\ \Rightarrow\ \delta\mathbf{X}^{\mathsf{T}}\
eq 0 \tag{16}\]  
leads to the following element formulation  
\[\begin{array}{l}\int\limits_{L}\int\limits_{A}\mathbf{N}^{\mathsf{T}}\mathbf{ \psi}^{\mathsf{T}}\mathbf{H}^{\mathsf{T}}\eta\mathbf{H}\Psi\frac{\partial \mathbf{X}}{\partial t^{2}}\mathbf{d}A\mathbf{d}\mathbf{z}+\int\limits_{L} \int\limits_{A}\mathbf{N}^{\mathsf{T}}\mathbf{\psi}^{\mathsf{T}}\mathbf{H}^{ \mathsf{T}}\mathbf{c}^{\mathsf{T}}\mathbf{E}\mathbf{c}\mathbf{H}\mathbf{\psi} \mathbf{d}A\mathbf{d}\mathbf{