## 2 Theory and Formulations
### Finite Element Modelling of the Axial FG Element with Material Damping  
Figure 2 shows a two-node finite element bar model \"\(e\)\", which has a total length of \(L\). The axial displacement is assumed in \(x\) direction along the bar, following the notation used in the literature for axial bar. The power law distribution determines the longitudinal variation of the material properties \(E(x)\), \(\rho(x)\) of the functionally graded bar. The function is defined as in Equation (1). \[P(x)=(P_{La}-P_{R})\Big{(}1-\frac{x}{L}\Big{)}^{na}+P_{R} \tag{1}\]  
where \(na\) is the power law exponent along the bar for material change. The axial displacement of any point of the bar in the \(x\) direction is  
\[u(x,z,t)=u_{0}(x,t) \tag{2}\]  
Considering the strain and stress relationships of the axial bar within the framework of [PERSON]'s law, the elastic strain energy of a finite element at any time is  
\[U=\frac{1}{2}\int\limits_{0}^{L}\int\limits_{A}E(x)\varepsilon_{xx}{}^{2}dAdx \tag{3}\]  
where \(\varepsilon_{xx}\) is the strain in \(x\) direction and for the area of the constant cross-section;  
\[A=\int_{A}dA \tag{4}\]  
The variation of [PERSON]'s modulus can be expressed as follows, obeying Equation (1):  
\[E(x)=(E_{L}-E_{R})\Big{(}1-\frac{x}{L}\Big{)}^{na}+E_{R} \tag{5}\]  
At any given moment, the axial displacement resulting in the kinetic energy of the Axially Functionally Graded bar in Cartesian coordinates is  
\[T=\frac{1}{2}\int\limits_{0}^{L}\int\limits_{A}\rho(x)\dot{u}^{2}dAdx \tag{6}\]  
where \(\dot{u}=\frac{du}{dt}\) is the velocity of the axial element. The variation of specific mass can be written as follows according to Equation (1):  
\[\rho(x)=(\rho_{L}-\rho_{R})\Big{(}1-\frac{x}{L}\Big{)}^{na}+\rho_{R} \tag{7}\]  
Figure 2: A two-node finite element bar model in local coordinate system. Each node has one degree of freedom for the bar element. To express the axial displacement \(u(x,t)\) of the element, the shape function and nodal displacements are used in matrix notation, as in Equation (8). \[\{u\}=[N]_{bx} \begin{Bmatrix}\dot{q}_{i}\\ q_{j}\end{Bmatrix}  \tag{8}\]  
The following extrapolation functions (Equation (9)) are used for the axial displacement along the bar:  
\[[N]_{bx}=\left\{1-\frac{x}{L}\ \ \frac{x}{L}\right\} \tag{9}\]  
By substituting Equation (8) into Equations (3) and (6), the energy functions can be rewritten as follows:  
\[U=\frac{1}{2}\left(\{q\}^{T}\Big{[}A_{bx}A\big{[}N^{\prime}\big{]}_{bx}^{T} \big{[}N^{\prime}\big{]}_{bx}\big{]}\{q\}\right) \tag{10}\]  
\[T=\frac{1}{2}\left(\{\dot{q}\}^{T}\left[\int\limits_{0}^{L}\rho(x)[N]_{bx}^{T} [N]_{bx}dx\right]\{\dot{q}\}\right) \tag{11}\]  
where  
\[\left[N^{\prime}\right]_{bx}=\frac{d[N]_{bx}}{dx} \tag{12}\]  
\[A_{bx}=\int\limits_{0}^{L}E(x)dx=\frac{(E_{L}-E_{R})L}{na+1}+E_{R}L \tag{13}\]  
Using the expressions between the displacement and velocity terms from Equations (10) and (11), the stiffness and mass matrix can be written as in Equation (14). \[ \begin{array}{l}[K_{A}]=\left[A_{bx}\ A[N^{\prime}]_{bx}^{T}[N^{\prime}]_{bx }\right]\\ [M_{A}]=\left[\int\limits_{0}^{L}\rho(x)[N]_{bx}^{T}[N]_{bx}dx\right]\end{array}  \tag{14}\]  
Hysteretic damping is a better choice for proportional damping if it is proportional to the stiffness matrix, rather than viscous damping, and it is only proportional to the stiffness matrix. However, hysteretic damping can also be proportional to the mass matrix. When the excitation is harmonic, the structural form of damping can only be obtained by replacing Young's modulus \(E\) with a complex one \(E(1+\eta)\). \(\eta\) is the material loss factor for isotropic homogenous bar. While deriving the element stiffness matrices, the complex matrix for material damping can be obtained as in Equation (15) [1,9]. \[ \begin{array}{l}[K]+i[H]\\ [H]=\eta[K]\end{array}  \tag{15}\]  
The axial damping coefficient is assumed to be proportional to Young modules by the loss factor.