## 2 Theory and Formulations
### Modelling Beam Element
By combining Equations (30) and (31) with Equations (24) and (27), the energy functions can be rewritten as follows:  
\[U=\frac{1}{2}\int_{0}^{L}\Bigl{[}A_{xx}\bigl{(}\left[N^{\prime}\right]_{bx} \cdot\{q\}_{u}\bigr{)}^{2}+D_{xx}(\left[N^{\prime\prime}\right]_{w}\cdot\{q\}_ {w})^{2}\Bigr{]}dx \tag{32}\]  
\[T=\frac{1}{2}\int_{0}^{L}\Bigl{[}I_{A}\bigl{(}\left[N\right]_{bx}\cdot\{\dot{q }\}_{u}\bigr{)}^{2}+I_{A}\bigl{(}\left[N\right]_{w}\cdot\{\dot{q}\}_{w}\bigr{)} ^{2}+I_{D}\bigl{(}\left[N^{\prime}\right]_{w}\cdot\{\dot{q}\}_{w}\bigr{)}^{2 }\Bigr{]}dx \tag{33}\]  
Equations (32) and (33) are expanded as  
\[U=\frac{1}{2}\biggl{[}\{q\}_{u}^{T}\biggl{(}A_{xx}\int_{0}^{L}\left[N^{\prime} \right]_{bx}\left[N^{\prime}\right]_{bx}dx\Bigr{)}\{q\}_{u}+\{q\}_{w}^{T} \biggl{(}D_{xx}\int_{0}^{L}\left[N^{\prime\prime}\right]_{w}\left[N^{\prime \prime}\right]_{w}dx\Bigr{)  
The local axis of AFG bar \"x\" lies along the z-axis of the beam and the global axis \(Z\) in the assemblage. The degrees of freedom at both the nodes of AFG bar and beam elements can therefore be transformed from local to global axes by means of the relation, as can be seen from Figure 4.  
The material change is defined along the bar by using Equation (1), as in Figures 2 and 3.  
In the assembly process, the degrees of freedom of beam are \(pb=3m+3\), \(m\) is the number of the finite elements and \(n=m+1\) is the number of the nodes. Considering the axial FG bar elements as support with the beam, total degrees of freedom of the system can be written as in Equation (37).  
\[p=3m+3+2 ma \tag{37}\]  
where last term \"_ma_\" belongs to the number of the elements of the axial FG bar and the degrees of freedom of axial FG bars are \(2 ma+2\). The number of elements of beam and each of axial FG bar is considered equal (64 elements). The assembly process of the elements of matrices can be derived in Equation (38) as follows:  
\[ \begin{array}{c}K(i,j)=K_{A}(i,j);i=1\ldots m,j=1\ldots m\\ K(m+i,\,m+j)=K_{B}(i,j),i=1\ldots pb,j=1\ldots pb\\ K(p+1-i,p+1-j)=K_{A}(i,j),i=1\ldots m,j=1\ldots m\end{array}  \tag{38}\]  
The assembly process of the intersection elements of matrices can be derived for the left and right supports, respectively, as in Equation (39), as follows:  
\[ \begin{array}{c}K(m,\,m+2)=K_{A}(m,\,m+1);\\ K(m+2,\,m)=K_{A}(m+1,\,m);\\ K(m+2,\,m+2)=K_{A}(m+1,\,m+1)+K(m+2,\,m+2);\\ K(4m+2,\,4m+4)=K_{A}(p,p\!-\!1);\\ K(4m+4,\,4m+2)=K_{A}(p\!-\!1,\,p);\\ K(4m+2,\,4m+2)=K_{A}(p,p)+K(4m+2,\,4m+2);\end{array}  \tag{39}\]  
where \(K\) is the global stiffness matrix, \(K_{A}\) is the axial FG bar stiffness matrix and \(K_{B}\) is the beam stiffness matrix. The same assembly processes is applied to mass and complex stiffness matrices of the system elements. After assembling the stiffness, complex stiffness and mass matrices, the global equations for a finite element model of the structure can be written by using the Lagrange equations. The set of differential equations governing the motion of the system is obtained considering the discrete (Equation (40)) and combined (Equation (41)) form of the system matrices as follows:  
\[[K_{B}]\{q_{B}\}+([K_{A}]+i[H])\{q_{A}\}+[M_{A}]\{\tilde{q}_{A}\}+[M_{B}]\{ \tilde{q}_{B}\}=\{F(t)\} \tag{40}\]  
Figure 4: Degrees of freedom used for beam and bar assembly in global coordinate system. \[[K]\{q\}+i[H]\{q\}+[M]\{\tilde{q}\}=\{F(t)\} \tag{41}\]  
The harmonic force acting at the centre of the beam is expressed as follows in Equation (42)  
\[\{F(t)\}=\{F_{\max}\}\cdot e^{i\omega t} \tag{42}\]  
\(\{F(t)\}\) is \(p\times 1\), a non-dimensional load vector. \(F_{\max}\) is the magnitude of force acting on the midpoint of the beam and is the only non-zero term in the force vector on the \(((p-1)/2)+1\) row. Time-dependent nodal displacements for steady-state vibrations of the system can be written as follows:  
\[\{q(t)\}=\{q_{o}\}e^{i\omega t} \tag{43}\]  
The elements of \(\{q_{o}\}\) are complex variables containing a phase angle.