Abeer M. Hessen*
| Mohammed Al-Shujairi | Ahmed F. Hamzah
© 2026 The authors. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).
OPEN ACCESS
This study examines how bi-directionally functionally graded (BDFG) beams respond to repeated moving harmonic loads. Unlike typical one-way functionally graded materials (FGMs), BDFG beams' material characteristics change along both the axial and thickness directions, thereby enhancing structural performance. This study aims to evaluate the free and forced vibration responses under various loading and support conditions. Using three beam theories—Euler–Bernoulli, Timoshenko, and higher-order shear deformation—the effects of material gradation, aspect ratio, boundary conditions, load velocity, excitation frequency, phase angle, and spacing between moving loads are investigated. The Lagrange equation was used to find the equations of motion, and the Newmark-Beta integration was used to solve them. The results indicate that the fundamental frequencies decrease with increasing gradient indices, reflecting the increased material compliance. Shear deformation has a significant influence on thick beams, and boundary conditions have a major influence on the dynamic response associated with load spacing and critical velocity. The stress patterns in BDFG beams are significantly different from those in homogeneous and one-way functionally graded beams. Moreover, dynamic deflections increase with load velocity and load spacing, reaching critical peaks at certain velocities.
bi-directionally functionally graded beam, moving harmonic loads, dynamic response, Shear deformation theory, Timoshenko beam theory, higher-order beam theory
Advanced composites known as functionally graded materials (FGMs) have qualities that change continuously in one or more directions, therefore removing the sharp edges seen in traditional laminated structures. Offering a great improvement over conventional composites, this smooth change lowers delamination and increases thermal and mechanical resistance [1, 2]. Natural systems, including teeth, bones, and bamboo, have inspired the creation of FGMs [3, 4] because they have structures that change from one place to another. Usually, ceramic-metal FGMs combine the strength of metals with the thermal resistance of ceramics, which makes them work better in places where it's hot or where there are chemicals [5]. Because of these benefits, FGMs are used a lot in things like cars, planes, nuclear power plants, and even in the medical field [6].
Most prior studies looked at one-dimensional FGMs, but more recently, people have started studying bi-directional functionally graded (BDFG) materials, where the properties change in both the axial and thickness directions. Because BDFGs offer better control over rigidity, vibration, and stress distribution, they are appealing for sophisticated structural uses. Several scholars have examined the dynamic and mechanical behavior of BDFG beams using several methods. For example, Simsek [7] used the Newmark-β technique to solve equations of motion for a 2D functionally graded Timoshenko beam under moving loads, whereas Deng and Cheng [8] looked at how gradients affected frequency and mode shapes. While Nejad et al. [9-11] used nonlocal elasticity theory to examine buckling, bending, and vibration of nanobeams, Wang et al. [12] studied natural frequencies of 2D functionally graded Euler–Bernoulli beams. While Shafiei and Kazemi [13] and Shafiei et al. [14] investigated buckling and vibration in porous micro- and nanobeams, Van Do et al. [15] examined static deflection and buckling of 2D functionally graded plates. Nguyen et al. [16] studied vibrational responses of functionally graded Timoshenko beams under moving loads While Trinh et al. [17] expanded the analysis to quasi-3D microbeams. Also studied have been nonlinear phenomena including bending, buckling, and free vibration of imperfect functionally graded nanobeams [18], rotating functionally graded beams [19], differential quadrature-based free vibration [20], and finite element techniques [21].
Notwithstanding these developments, there are still a lot of problems. It is still hard to precisely identify mode shapes of BDFG beams, hence many research relies on approximations employing homogeneous modes [22]. Moreover, most earlier studies have concentrated on free vibration or reduced load situations. Forced vibration of BDFG beams exposed to repeated moving harmonic loads has not received much attention, especially with different boundary conditions [23].
The current research aims to fill this vacuum by examining the free and forced vibration responses of BDFG beams exposed to two sequential moving harmonic loads. Three beam theories, Euler–Bernoulli, Timoshenko, and higher-order shear deformation, are used to record shear effects and geometric influences. Lagrange's formulation is used to get the equations of motion, and the Newmark-β integration method is used to solve them. There are three main things that this work does:
(i) It gives a comparison of different beam theories for BDFG structures.
(ii) It looks at how things like gradient indices, aspect ratio, boundary conditions, load velocity, and excitation frequency affect dynamic behavior.
(iii) It gives some ideas on how to make BDFG beams better for advanced engineering uses.
2.1 Displacement field and assumptions
The BDFG beams experienced repeated moving harmonic loads. The displacement field, assumptions, material gradation, energy equations, approximation techniques, and last motion equations are all part of the development. Two concentrated loads (P1, P2) are put on the BDFG beam, which has a length (L), a width (b), and a thickness (h) as shown in Figure 1. L, b, and h are chosen to match the length, breadth, and thickness of the BDFG beam. Along with the fundamental coordinates (x,y,z), and under the premise that the BDFG beam is fixed in four boundary conditions: the clamped–clamped (CC), the clamped–simply supported (CS), the simply supported (SS), and the clamped-free (CF). Along its longitudinal (x) and thickness (z) directions, the BDFG beam, a composite of metals and ceramics, has constantly changing effective material properties: Young's modulus (E), Poisson's ratio (ν), shear modulus (G), and mass density (ρ).
Figure 1. The geometry of a beam made of a bi-directionally functionally graded beam subjected to two successive moving harmonic loads
2.2 Material gradation
The BDFG beam is made of a mix of ceramic and metal. Its properties change smoothly along the x-axis and thickness (z-axis) directions. Young's modulus E, Poisson's ratio υ, shear modulus G, and mass density ρ are the features of this theory that make use of the same Eqs. (1)-(3):
$P(x, z)=P_{L B} e^{k_1 \alpha(x)+k_3 \beta(z)}$ (1)
$\alpha(x)=\frac{x}{L}+\frac{1}{2}$ (2)
$\beta(z)=\frac{z}{h}+\frac{1}{2}$ (3)
Here, P(x,z) represents the material properties (e.g., E, ν, ρ, G). PLB denotes the reference material property evaluated at the coordinates (L/2; h/2), as illustrated in Figure 2. k1 and k3 are the gradient indices that determine the variation pattern of the material along the x and z axes, respectively. Setting k1 and k3 to zero makes the beam homogeneous, characterized by a constant modulus of elasticity, denoted as E = ELB.
Figure 2. The variation of the elastic modulus of a bi-directionally functionally graded beam
2.3 Energy expressions
The potential energy of the beam is derived from strain energy due to bending, shear deformation, and external moving harmonic loads. The kinetic energy is formulated in terms of the velocity components of the displacement field. The displacement field in higher-order shear deformation theory (HOSDT) is:
$u_x(x, z, t) 1=u_0(x, t)-z \frac{\partial \omega}{\partial x}+A(z) \gamma(x, t)$ (4)
$u_y(x, z, t) 1=0$ (5)
$u_z(x, z, t) 1=\omega(x, t)$ (6)
where, u0 and ω are the axial and transverse displacements of any point on the neutral axis.
A(z) is a function of the stress distribution along the thickness of the beam that characterizes the transverse shear. The displacement fields of different beam theories can be obtained by choosing the displacement field. A(z) as follows:
• Euler-Bernoulli beam theory (EBBT) or classical beam theory: A(z) = 0
• First order shear deformation beam theory (FSDBT) or Timoshenko beam theory (TBT): A(z) = z
• Parabolic shear deformation beam theory (PSDBT) (Reddy): $A(z)=z\left(1-\frac{4 z^2}{3 h^2}\right)$
$\gamma(x, t) 1=\frac{\partial \omega}{\partial x}-\phi(x, t)$ (7)
The potential energy equation for parabolic shear deformation beam theory.
$U_{i n}=\frac{1}{2} \int_V\left(\sigma_x \varepsilon_x+\tau_{x z} \gamma_{x z}\right) d V$ (8)
The equation after simplification becomes:
$U_{i n}=\frac{1}{2} \int_0^L \int_A\left(E(x, z) 1 \varepsilon_x^2+G(x, z) \gamma_{x z}^2\right) d A d x$ (9)
Hint:
$\begin{gathered}A_{x x}=\int_A E_{L B} e^{k_3 \beta(z)}(1) d A \rightarrow B_{x x}=\int_A E_{L B} e^{k_3 \beta(z)}(z) d A \\ D_{x x}=\int_A E_{L B} e^{k_3 \beta(z)}\left(z^2\right) d A \rightarrow E_{x x}=\int_A A(z) E_{L B} e^{k_3 \beta(z)}(1) d A \\ F_{x x}=\int_A A(z) 1 E_{L B} e^{k_3 \beta(z)}(z) d A \rightarrow \\ H_{x x}=\int_A E_{L B} e^{k_3 \beta(z)}(A(z))^2 d A 1 B_{x z}=\int_A G(x, z)\left(A^{\prime}(z)\right)^2 d A\end{gathered}$
The equation after simplification becomes:
$U_{i n}=\left[\begin{array}{l}\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} A_{x x}\left(u^{\prime}\right)^2 d x \\ +\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} D_{x x}\left(\omega^{\prime \prime}\right)^2 d x \\ +\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} H_{x x}\left(\gamma^{\prime}\right)^2 d x \\ -\int_0^L e^{k_1 \alpha(x)} B_{x x} u^{\prime} \omega^{\prime \prime} d x \\ +\int_0^L e^{k_1 \alpha(x)} E_{x x} u^{\prime} \gamma^{\prime} d x \\ -\int_0^L e^{k_1 \alpha(x)} F_{x x} \omega^{\prime \prime} \gamma^{\prime} d x \\ +\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} B_{x z} \gamma^2 d x\end{array}\right]$ (10)
The kinetic energy of the BDFG beam can be written as:
$K_e=\frac{1}{2} \int_V \rho(x, z)\left[\dot{u}_x^2+\dot{u}_y^2+\dot{u}_z^2\right] d V$ (11)
The kinetic energy of the BDFG beam can be written as:
$K=\left[\begin{array}{c}\frac{1}{2} \int_0^L I_A(\dot{u})^2 d x+\frac{1}{2} \int_0^L I_D\left(\dot{\omega}^{\prime}\right)^2 d x+\frac{1}{2} \int_0^L I_H(\dot{\gamma})^2 d x \\ -\int_0^L I_B(\dot{u})\left(\dot{\omega}^{\prime}\right) d x \\ +\int_0^L I_E(\dot{u})(\dot{\gamma}) d x-\int_0^L I_F\left(\dot{\omega}^{\prime}\right)(\dot{\gamma}) d x \\ +\frac{1}{2} \int_0^L I_A \dot{\omega}^2 d x\end{array}\right]$ (12)
The inertia coefficients appearing in Eq. (12) are defined as:
$\begin{gathered}\left(I_A, I_B, I_D\right)=\int_A \rho(x, z)\left(1, z, z^2\right) d A \\ \left(I_E, I_F, I_H\right)=\int_A \rho(x, z) A(z)\left(1, z, z^2\right) d A\end{gathered}$
2.4 Galerkin approximation
To apply Lagrange's equations, the unknown displacement functions w0(x,t), u0(x,t), and γ(x,t) can be expressed in terms of generalized time-dependent coordinates, and the admissible functions must satisfy basic geometric boundary conditions. Hence, the displacement functions can be expressed as follows:
$\omega(x, t)=\sum_{m=1}^N A_m(t) \phi_{m w}(x), \phi_{m w}(x)=g_w(x) x^{m-1}$ (13)
$u(x, t)=\sum_{m=1}^N B_m(t) \alpha_{m u}(x), \alpha_{m u}(x)=g_u(x) x^{m-1}$ (14)
$\gamma(x, t)=\sum_{m=1}^N C_m(t) \beta_{m \gamma}(x), \beta_{m \gamma}(x)=g_\gamma(x) x^{m-1}$ (15)
where, Am(t), Bm(t), Cm(t) are unknown generalized coordinates to be determined. gw(x), gu(x), gγ(x) are complementary functions that meet the fundamental boundary constraints.
$g_w(x)=\left(x+\frac{L}{2}\right)^{P_w}\left(x-\frac{L}{2}\right)^{q_w}$ (16)
$g_u(x)=\left(x+\frac{L}{2}\right)^{P_u}\left(x-\frac{L}{2}\right)^{q_u}$ (17)
$g_\gamma(x)=\left(x+\frac{L}{2}\right)^{P_\phi}\left(x-\frac{L}{2}\right)^{q_\phi}$ (18)
where, $P_i$ and $q_i(i=u, w, \phi)$ are the boundary exponents that are related to the boundary conditions. The related values of $P_i$ and $q_i$ are provided in Table for the 11 considered boundary conditions.
Table 1 defines the various boundary conditions used for higher-order beam theory and the corresponding values of the boundary-related exponents. For the SS, the values Pw = 1, Pu = 1, and Pφ = 0 are assigned at both ends, indicating that the transverse displacement and rotation are restricted while axial displacement conditions are relaxed. In the CC case, the values are Pw = 2, Pu = 1, and Pφ = 1 at both ends, with all corresponding q-values equal to (qw = 2, qu = 1, qφ = 1), representing complete restrictions on axial, transverse, and rotational displacements. For the CS, the left end adopts clamped conditions (Pw = 2, Pu = 1, Pφ = 1) while the right end corresponds to SS conditions (qw = 1, qu = 0, qφ = 0), reflecting the mixed boundary effect. Finally, for the CF, the left end is clamped with full restrictions (Pw = 2, Pu = 1, Pφ = 1), while the right end is free (qu = qw =qφ = 0), meaning no constraints are applied at that end. This classification provides the foundation for constructing the admissible functions used in the Galerkin approximation to enforce correct boundary behavior in the beam vibration analysis.
Table 1. Various values of the boundary exponents for different beam types
|
Boundary Conditions of Higher-Order Beam Theory |
Left End |
Right End |
||||
|
Pw |
Pu |
Pφ |
qw |
qu |
qφ |
|
|
Simply supported beam |
1 |
1 |
0 |
1 |
0 |
0 |
|
Clamped–clamped beam |
2 |
1 |
1 |
2 |
1 |
1 |
|
Clamped–simply supported beam |
2 |
1 |
1 |
1 |
0 |
0 |
|
Clamped-free beam |
2 |
1 |
1 |
0 |
0 |
0 |
$U_{i n t}=\left[\begin{array}{c}\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} A_{x x}\left(B_m \alpha_{m u}^{\prime}\right)^2 d x+\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} D_{x x}\left(A_m \phi_{m w}^{\prime \prime}\right)^2 d x \\ +\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} H_{x x}\left(C_m \beta_{m \gamma}^{\prime}\right)^2 d x-\int_0^L e^{k_1 \alpha(x)} B_{x x}\left(B_m \alpha_{m u}^{\prime}\right)\left(A_m \phi_{m w}^{\prime \prime}\right) d x \\ +\int_0^L e^{k_1 \alpha(x)} E_{x x}\left(B_m \alpha_m^{\prime}\right)\left(C_m \beta_{m \gamma}^{\prime}\right) d x-\int_0^L e^{k_1 \alpha(x)} F_{x x}\left(A_m \phi_{m w}^{\prime \prime}\right)\left(C_m \beta_{m \gamma}^{\prime}\right) d x \\ +\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} B_{x z}\left(C_m \beta_{m \gamma}\right)^2 d x\end{array}\right]$ (13)
$K e=\left[\begin{array}{c}\frac{1}{2} \int_0^L I_A\left(\dot{B}_m \cdot \alpha_{m u}\right)^2 d x+\frac{1}{2} \int_0^L I_D\left(\dot{A}_m \cdot \phi_{m w}^{\prime}\right)^2 d x+\frac{1}{2} \int_0^L I_H\left(\dot{C}_m \cdot \beta_{m \gamma}\right) d x \\ -\int_0^L I_B\left(\dot{B}_m \cdot \alpha_{m u}\right)\left(\dot{A}_m \phi_{m w}^{\prime}\right) d x+\int_0^L I_E\left(\dot{B}_m \cdot \alpha_{m u}\right)\left(\dot{C}_m \beta_{m \gamma}\right) d x \\ -\int_0^L I_F\left(\dot{A}_m \cdot \phi_m^{\prime}\right)\left(\dot{C}_m \beta_{m \gamma}\right) d x+\frac{1}{2} \int_0^L I_A\left(\dot{A}_m \cdot \phi_{m w}\right)^2 d x\end{array}\right]$ (14)
Using the Lagrange Method:
$L=\left[\begin{array}{l}\frac{1}{2} \int_0^L I_A\left(\dot{B}_m \cdot \alpha_{m u}\right)^2 d x+\frac{1}{2} \int_0^L I_D\left(\dot{A}_m \cdot \phi_{m v}^{\prime}\right)^2 d x+\frac{1}{2} \int_0^L I_H\left(\dot{C}_m \cdot \beta_{m \gamma}\right)^2 d x \\ -\int_0^L I_B\left(\dot{B}_m \cdot \alpha_{m u}\right)\left(\dot{A}_m \phi_{m v}^{\prime}\right) d x+\int_0^L I_E\left(\dot{B}_m \cdot \alpha_{m u}\right)\left(\dot{C}_m \beta_{m \gamma}\right) d x-\int_0^L I_F\left(\dot{A}_m \cdot \phi_m^{\prime}\right)\left(\dot{C}_m \beta_{m \gamma}\right) d x \\ +\frac{1}{2} \int_0^L I_A\left(\dot{A}_m \cdot \phi_{m v}\right)^2 d x-\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} A_{x x}\left(B_m \alpha_{m u}^{\prime}\right)^2 d x-\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} D_{x x}\left(A_m \phi_{m v}^{\prime \prime}\right)^2 d x \\ -\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} H_{x x}\left(C_m \beta_{m \gamma}^{\prime}\right)^2 d x+\int_0^L e^{k_1 \alpha(x)} B_{x x}\left(B_m \alpha_{m u}^{\prime}\right)\left(A_m \phi_{m w}^{\prime \prime}\right) d x-\int_0^L e^{k_1 \alpha(x)} E_{x x}\left(B_m \alpha_m^{\prime}\right)\left(C_m \beta_{m \gamma}^{\prime}\right) d x \\ +\int_0^L e^{k_1 \alpha(x)} F_{x x}\left(A_m \phi_{m w}^{\prime \prime}\right)\left(C_m \beta_{m \gamma}^{\prime}\right) d x-\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} B_{x z}\left(C_m \beta_{m \gamma}\right)^2 d x\end{array}\right]$ (15)
$\begin{aligned} & \frac{\partial L}{\partial A_k}-\frac{d}{d t}\left(\frac{\partial L}{\partial \dot{A}_k}\right) \\ & =0-\frac{1}{2} \int_0^L e^{k_1 \alpha(x)} D_{x x} . 2\left(A_m \phi_{m w}^{\prime \prime}\right) \phi_k^{\prime \prime} d x \\ & +\int_0^L e^{k_1 \alpha(x)} B_{x x}\left(B_m \alpha_{m u}^{\prime}\right) \phi_k^{\prime \prime} d x \\ & +\int_0^L e^{k_1 \alpha(x)} F_{x x}\left(C_m \beta_{m \gamma}^{\prime}\right) \phi_k^{\prime \prime} d x \\ & -\frac{d}{d t}\left[\frac{1}{2} \int_0^L I_D . 2\left(\dot{A}_m . \phi_{m w}^{\prime}\right) \phi_k^{\prime} d x\right. \\ & -\int_0^L I_B\left(\dot{B}_m . \alpha_{m u}\right) \phi_k^{\prime} d x-\int_0^L I_F\left(\dot{C}_m \beta_{m \gamma}\right) \phi_k^{\prime} d x \\ & \left.+\frac{1}{2} \int_0^L I_A . 2\left(\dot{A}_m \phi_m\right) \phi_k d x\right]=0\end{aligned}$ (16)
The equation after simplification becomes:
$\begin{aligned} & -\int_0^L e^{k_1 \alpha(x)} D_{x x} \phi_{m w}^{\prime \prime} \phi_k^{\prime \prime} A_m d x \\ & +\int_0^L e^{k_1 \alpha(x)} B_{x x} \alpha_{m u}^{\prime} \phi_k^{\prime \prime} B_m d x \\ & +\int_0^L e^{k_1 \alpha(x)} F_{x x} \beta_{m \gamma}^{\prime} \phi_k^{\prime \prime} C_m d x \\ & -\int_0^L I_D \phi_{m w}^{\prime} \phi_k^{\prime} \ddot{A}_m d x \\ & +\int_0^L I_B \alpha_{m u} \phi_k^{\prime} \ddot{B}_m d x \\ & +\int_0^L I_F \beta_{m \gamma} \phi_k^{\prime} \ddot{C}_m d x \\ & -\int_0^L I_A \phi_m \phi_k \ddot{A}_m d x=0\end{aligned}$ (17)
We get the terms of stiffness and mass matrix
$\begin{aligned} & K_1=-\int_0^L e^{k_1 \alpha(x)} D_{x x} \phi_{m w}^{\prime \prime} \phi_k^{\prime \prime} d x \\ & K_2=\int_0^L e^{k_1 \alpha(x)} B_{x x} \alpha_{m u}^{\prime \prime \prime} d x \\ & K_3=\int_0^L e^{k_1 \alpha(x)} F_{x x} \beta_{m \gamma}^{\prime} \phi^{\prime \prime} d x\end{aligned}$
$\begin{gathered}M_1=-\int_0^L I_D \phi_{m w}^{\prime} \phi_k^{\prime} d x-\int_0^L I_A \phi_m \phi_k d x \\ M_2=\int_0^L I_B \alpha_{m u} \phi_k^{\prime} d x \\ M_3=\int_0^L I_F \beta_{m \gamma} \phi_k^{\prime} d x\end{gathered}$
$\begin{aligned} & \frac{\partial L}{\partial B_k}-\frac{d}{d t}\left(\frac{\partial L}{\partial \dot{B}_k}\right)=0 \\ & -\int_0^L e^{k_1 \alpha(x)} A_{x x}\left(B_m \alpha_{m u}^{\prime}\right) \alpha_k^{\prime} d x \\ & +\int_0^L e^{k_1 \alpha(x)} B_{x x}\left(A_m \phi_{m w}^{\prime \prime}\right) \alpha_k^{\prime} d x \\ & -\int_0^L e^{k_1 \alpha(x)} E_{x x}\left(C_m \beta_{m \gamma}^{\prime}\right) \alpha_k^{\prime} d x \\ & -\int_0^L I_A\left(\ddot{B}_m \cdot \alpha_{m u}\right) \alpha_k d x \\ & +\int_0^L I_B\left(\ddot{A}_m \phi_{m w}^{\prime}\right) \alpha_k d x \\ & -\int_0^L I_E\left(\ddot{C}_m \beta_{m \gamma}\right) \alpha_k d x=0\end{aligned}$ (18)
From Eq. (17) we get that:
$\begin{gathered}K_4=\int_0^L e^{k_1 \alpha(x)} B_{x x} \phi_{m w}^{\prime \prime} \alpha_k^{\prime} d x \\ K_5=-\int_0^L e^{k_1 \alpha(x)} A_{x x} \alpha_{m u}^{\prime} \alpha_k^{\prime} d x \\ K_6=-\int_0^L e^{k_1 \alpha(x)} E_{x x} \beta_{m \gamma}^{\prime} \alpha_k^{\prime} d x \\ M_4=\int_0^L I_B \phi_{m w}^{\prime} \alpha_k d x \\ M_5=-\int_0^L I_A \alpha_{m u} \alpha_k d x \\ M_6=-\int_0^L I_E \beta_{m \gamma} \alpha_k d x\end{gathered}$
$\begin{aligned} & \frac{\partial L}{\partial C_k}-\frac{d}{d t}\left(\frac{\partial L}{\partial \dot{C}_k}\right)=0 \\ & -\int_0^L e^{k_1 \alpha(x)} H_{x x}\left(C_m \beta_{m \gamma}^{\prime}\right) \beta_k^{\prime} d x \\ & -\int_0^L e^{k_1 \alpha(x)} E_{x x}\left(B_m \alpha_m^{\prime}\right) \beta_k^{\prime} d x \\ & +\int_0^L e^{k_1 \alpha(x)} F_{x x}\left(A_m \phi_{m w}^{\prime \prime}\right) \beta_k^{\prime} d x \\ & -\int_0^L e^{k_1 \alpha(x)} B_{x z}\left(C_m \beta_{m \gamma}\right) \beta_k d x \\ & -\int_0^L I_H\left(\ddot{C}_m \cdot \beta_{m \gamma}\right) \beta_k d x \\ & -\int_0^L I_E\left(\ddot{B}_m \cdot \alpha_{m u}\right) \beta_k d x \\ & +\int_0^L I_F\left(\ddot{A}_m \cdot \phi_m^{\prime}\right) \beta_k d x=0\end{aligned}$ (19)
From Eq. (18) we get that:
$\begin{gathered}K_7=\int_0^L e^{k_1 \alpha(x)} F_{x x} \phi_{m w}^{\prime \prime} \beta_k^{\prime} d x \\ K_8=-\int_0^L e^{k_1 \alpha(x)} E_{x x} \alpha_m^{\prime} \beta_k^{\prime} d x \\ K_9=-\int_0^L e^{k_1 \alpha(x)} H_{x x} \beta_{m \gamma}^{\prime} \beta_k^{\prime} d x-\int_0^L e^{k_1 \alpha(x)} B_{x z} \beta_{m \gamma} \beta_k d x \\ M_7=\int_0^L I_F \phi_m^{\prime} \beta_k d x \\ M_8=-\int_0^L I_E \alpha_{m u} \beta_k d x \\ M_9=-\int_0^L I_H \beta_{m \gamma} \beta_k d x\end{gathered}$
$\left[\begin{array}{lll}K_1 & K_2 & K_3 \\ K_4 & K_5 & K_6 \\ K_7 & K_8 & K_9\end{array}\right]\left\{\begin{array}{l}A_m \\ B_m \\ C_m\end{array}\right\}+\left[\begin{array}{lll}M_1 & M_2 & M_3 \\ M_4 & M_5 & M_6 \\ M_7 & M_8 & M_9\end{array}\right]\left\{\begin{array}{l}\ddot{A}_m \\ \ddot{B}_m \\ \ddot{C}_m\end{array}\right\}=\left\{\begin{array}{l}0 \\ 0 \\ 0\end{array}\right\}$ (20)
where, $K_1-K_7$ present stiffness matrices, $M_1-M_7$ present mass matrices. For free vibration analysis, the time-dependent generalized coordinates can be written as:
$A_m(t)=\bar{A}_m e^{i \omega t}, B_m(t)=\bar{B}_m e^{i \omega t}, C_m(t)=\bar{C}_m e^{i \omega t}$
where, $i=\sqrt{-1}$ and $\omega$ is the natural frequency. Setting the external load vector to zero, and substituting equations above into Eq. (20) yields the following frequency equation:
$\left[\begin{array}{lll}K_1 & K_2 & K_3 \\ K_4 & K_5 & K_6 \\ K_7 & K_8 & K_9\end{array}\right]-\omega\left[\begin{array}{lll}M_1 & M_2 & M_3 \\ M_4 & M_5 & M_6 \\ M_7 & M_8 & M_9\end{array}\right]\left\{\begin{array}{l}\bar{A} \\ \bar{B} \\ \bar{C}\end{array}\right\}=\left\{\begin{array}{l}0 \\ 0 \\ 0\end{array}\right\}$ (21)
The natural frequencies can be obtained from the condition that the determinant of the coefficient matrix is zero.
2.5 Formulation of forced vibration
Two consecutive moving harmonic loads are used to examine the BDFG beams' forced vibration. Lagrange's equations are used to get the governing equations of motion in generalized coordinates after the displacement field is plugged into the energy expressions. Load velocity and excitation frequency affect the time-dependent harmonic function that describes the external excitation.
The implicit Newmark-β method is used to numerically solve the ensuing system of coupled equations, therefore producing stable and precise time responses. This expression helps to evaluate how load factors, boundary conditions, and gradient indices influence the dynamic response of BDFG beams. The following ideas underpin this research: (i) The beam is first motionless; hence, the initial conditions are zero; (ii) The moving load's velocity is constant, and it stays in touch with the beam during the excitation; and (iii) The moving load's inertial effects are deemed unimportant.
2.6 Equations of motion
Applying Lagrange's equations called for finding the potential energy of two nearby moving harmonic loads moving at the same speed. The possibilities of the outside loads were then stated as follows:
$U_{ {ext }}=\left\{\begin{array}{cc}-P_1(t) w\left(x_{p 1}, t\right) & { if } 0 \leq t \leq t_1=\frac{d}{v_p} \\ -P_1(t) w\left(x_{p 1}, t\right)-P_2(t) w\left(x_{p 2}, t\right) & \text { if } t_1 \leq t \leq t_2=\frac{L}{v_p} \\ -P_2(t) w\left(x_{p 2}, t\right) & \text { if } t_2 \leq t \leq t_1+t_2 \\ 0 & f t_1+t_2 \leq t\end{array}\right\}$ (22)
where,
$P_1(t)=P_0 \sin \left(\Omega_1 t+\phi_1\right)$ (23)
$P_2(t)=P_0 \sin \left(\Omega_2 t+\phi_2\right)$ (24)
The equation for moving harmonic loads is as follows: $P_0$ is the amplitude, $\Omega_1$ and $\Omega_2$ are the excitation frequencies, $\phi_1$ and $\phi_2$ are the phase angles, $t_1$ is the time when the load $P_1(t)$ first enters the beam, $t_2$ is the time when the load $P_2(t)$ leaves the beam, and $x_{p 1} x_{p 1}$ and $x_{p 2}$ are the positions of the first and second moving harmonic loads at any given instant.
$x_{p 1}=v_p t-L / 2$ (25)
$x_{p 2}=v_p t-L / 2-d$ (26)
where, $d$ represents the distance between two consecutive moving loads, and $v_p$ denotes the velocity of the moving harmonic loads in the axial direction. To apply Lagrange's equations, the unknown displacement functions $w_0(x, t)$, $u_0(x, t)$ and $\phi(x, t)$ can be expressed in terms of timedependent generalized coordinates and admissible functions that must adhere to the essential (geometric) boundary conditions. Consequently, the displacement functions may be expressed as the following series.
$\begin{aligned} & \omega(x, t)=\sum_{m=1}^N A_m(t) \phi_{m w}(x), \\ & \phi_{m w}(x)=g_w(x) x^{m-1}\end{aligned}$ (27)
$\begin{aligned} & u_0(x, t)=\sum_{m=1}^N B_m(t) \alpha_{m u}(x), \\ & \alpha_{m u}(x)=g_u(x) x^{m-1}\end{aligned}$ (28)
$\begin{aligned} & \phi(x, t)=\sum_{m=1}^N C_m(t) \beta_{m \phi}(x), \\ & \beta_{m \phi}(x)=g_\phi(x) x^{m-1}\end{aligned}$ (29)
where, $A_m(t), B_m(t), C_m(t)$ are known generalized coordinates to be determined, and $g_{u(x)}, g_{u(x)}, g_{d(x)}$ are auxiliary functions used to satisfy essential boundary conditions. The auxiliary functions are written in the following form:
$\begin{aligned} g_w(x) & =\left(x+\frac{L}{2}\right)^{P_w}\left(x-\frac{L}{2}\right)^{q_w} \\ g_u(x) & =\left(x+\frac{L}{2}\right)^{P_u}\left(x-\frac{L}{2}\right)^{q_u} \\ g_\phi(x) & =\left(x+\frac{L}{2}\right)^{P_\phi}\left(x-\frac{L}{2}\right)^{q_\phi}\end{aligned}$
where, $P_i$ and $q_i(i=u, w, \phi)$ are the boundary exponents, which are contingent on the boundary conditions, with the corresponding $P_i$ and $q_i$ values provided in Table 1. As previously mentioned, the equations of motion for this problem are derived using Lagrange's equations:
$\frac{d}{d t}\left(\frac{\partial K_e}{\partial \dot{A}_k}\right)+\frac{\partial U_{i n t}}{\partial A_k \frac{\partial U_{e x t}}{\partial A_k}}, k=1,2,3, \ldots, 3 N$ (30)
where, $A_k=A_k, k=1,2,3, \ldots \ldots \ldots, N, A_k=B_{k-N}, k= N+1, \ldots \ldots ., 2 N, A_k=C_{k-2 N}, k=2 N+1, \ldots \ldots, 3 N$.
Applying the Lagrange equations yields the following couple of equations of motion:
$\begin{aligned} & {\left[\begin{array}{ccc}K_1 & 0 & K_3 \\ 0 & K_5 & K_6 \\ K_7 & K_8 & K_9\end{array}\right]\left\{\begin{array}{l}A(t) \\ B(t) \\ C(t)\end{array}\right\}} +\left[\begin{array}{ccc}M_1 & 0 & 0 \\ 0 & M_5 & M_6 \\ 0 & M_8 & M_9\end{array}\right]\left\{\begin{array}{l}\ddot{A}(t) \\ \ddot{B}(t) \\ \ddot{C}(t)\end{array}\right\}=\left\{\begin{array}{c}F(t) \\ 0 \\ 0\end{array}\right\}\end{aligned}$ (31)
$K_1-K_9$ represent the stiffness matrices, $M_1-M_9, K_i$ and $M_i$ denote the mass matrices, and $F(t)$ represents the timedependent generalized load vector from two moving harmonic loads. The implicit Newmark- $\beta$ time integration method can solve the equations of motion. It should be noted here that the mass and stiffness matrices are symmetric and in size $N \times N$. Eq. (31) can be changed into a short matrix form, which looks like:
$[K]\{A(t)\}+[M]\{\ddot{A}(t)\}=\{Q(t)\}$ (32)
where, $\{A(t)\}=\{\{A(t)\},\{B(t)\},\{C(t)\}\}^T$ and $\{Q(t)\}= \{\{F(t)\},\{0\},\{0\}\}^T$.
Then, the Newmark time integration method is used to solve the time-dependent multiple equations. Finally, the beam's displacements, speeds, and accelerations at the time and point in question are found for any time $t$ between $0 \leq t \leq L / v_p$.
In Eq. (32), the following abbreviations have been introduced:
$\begin{gathered}K_1=-A_{x z} \int_{-L / 2}^{L / 2}\left(\theta_m\right)^{\prime}\left(\theta_k\right)^{\prime} d x \\ K_2=-A_{x z} \int_{-L / 2}^{L / 2}\left(\beta_m\right)\left(\theta_k\right)^{\prime} d x \\ K_3=-A_{x z} \int_{-L / 2}^{L / 2}\left(\alpha_m\right)^{\prime}\left(\alpha_k\right)^{\prime} d x \\ K_4=-B_{x x} \int_{-L / 2}^{L / 2}\left(\beta_m\right)^{\prime}\left(\alpha_k\right)^{\prime} d x \\ K_5=A_{x z} \int_{-L / 2}^{L / 2}\left(\theta_m\right)^{\prime}\left(\beta_k\right) d x \\ K_6=B_{x x} \int_{-L / 2}^{L / 2}\left(\alpha_m\right)^{\prime}\left(\beta_k\right)^{\prime} d x \\ M_1=-I_A \int_{-L / 2}^{L / 2} \theta_m \theta_k d x \\ M_2=-I_A \int_{-L / 2}^{L / 2} \alpha_m \alpha_k d x \\ M_3=-I_B \int_{-L / 2}^{L / 2} \beta_m \alpha_k d x \\ M_4=-I_B \int_{-L / 2}^{L / 2} \beta_m \beta_k d x \\ M_5=-I_D \int_{-L / 2}^{L / 2} \beta_m \beta_k d x \\ F=-P_1(t) \theta_k\left(x_{p 1}\right), \text { if } 0 \leq t \leq t_1=\frac{d}{v_p} \\ F=-P_1(t) \theta_k\left(x_{p 1}\right)-P_2(t) \theta_k\left(x_{p 2}\right), \text { if } t_1 \leq t \leq t_2=\frac{L}{v_p} \\ F=-P_2(t) \theta_k\left(x_{p 2}\right), \text { if } t_2 \leq t \leq t_1+t_2\end{gathered}$
where, $m, k=1,2, \ldots \ldots, N$.
Using numerical analysis, we investigate here the free and forced vibrations of BDFG beams. We investigate how the time history, excitation frequencies of the moving harmonic loads, and boundary conditions on these vibration responses are influenced by material composition (or gradient indices), moving load velocity (Vp), and the separation between two coordinate moving loads.
The beam's geometrical and material characteristics are as follows: $E_{L B}=210\mathrm{Gp}$, $v_{L B}=0.3$, $\rho_{L B}=7850\mathrm{~kg} / \mathrm{m}^3$, $b=0.5 \mathrm{~m}$, h = 0.1 and the length of the beam (L) is changed to see what happens with the shear deformation or (Aspect ratio, $A.R=L / h$). The shear correction factor is assumed to have the constant value of $K_s=5 / 6$ for the rectangular cross-section. The analysis considers SS, CC, CS, and CF boundary conditions, where S, C, and F represent simple, clamped, and free edges, respectively. The exponents $P_w=q_w=2$ should be selected for CC end conditions for the Euler-Bernoulli beam, and the boundary exponent $P_\phi$ is not applicable for EBBT. respect to the appropriate boundary exponents provided in Table 1. The dynamic deflections in the results of forced vibration are normalized by the static deflection of the SS homogeneous beam with $E=E_{L B}$ as:
$w_{s t}=\frac{P . L^3}{48 E_{L B} I}$ (33)
The dimensionless time $t^*$ is defined by
$t^*=\frac{x_p}{L}=\frac{v_p t}{L}-\frac{1}{2}$
From the equation above, it can be seen that the moving load is on the left side of the beam when $t^*$ = −0.5 and on the right side of the beam when $t^*$ = 0.5.
3.1 Responses to free vibration
Table 2 gives the dimensionless fundamental frequencies (ω) of BDFG beams with SS end conditions for two distinct aspect ratios (L/h = 5, 20) and various material gradient indices (k1, k3). The analyses are conducted for EBBT and TBT to examine the shear effect.
In the initial phase of validation studies, the dimensionless fundamental frequencies from the current formulation are contrasted [24]. In Table 2, the dimensionless fundamental frequencies of the functionally sorted bidirectional beam (SS) were calculated using (EBBT and TBT) and compared with the results of the study by Kapuria et al. [24] for slenderness ratio (i.e., L/h = 5, 20), and gradient index values (i.e., k1, k3 = 0, 0.2, 0.4, 0.6, 0.8, 1). For the comparison, the following frequency parameter in dimensionless form is used:
$\bar{\omega}=\frac{\omega L^2}{h} \sqrt{\frac{P_{L B}}{E_{L B}}}$ (34)
The current frequency values align with the prior findings in the study by Kapuria et al. [24]. As shown in Tables 2 and 3, TBT yields lower natural frequencies than the EBBT or HOSDT because TBT accounts for shear effects. Shear deformation reduces the beam's stiffness. The effect of shear deformation is significant for thick beams (L/h ≤ 10) but negligible for thin beams (L/h ≥ 20), leading to nearly identical frequency results for TBT and EBBT in the latter case.
Table 2. Dimensionless frequency of simply supported (SS) bi-directional functionally graded materials in free vibration
|
Beam Type |
Gradient Index, k1 |
Beam Theory |
L/h = 5 |
L/h = 20 |
||||||||||
|
Gradient Index, k3 |
Gradient Index, k3 |
|||||||||||||
|
k3 = 0 |
k3 = 0.2 |
k3 = 0.4 |
k3 = 0.6 |
k3 = 0.8 |
k3 = 1 |
k3 = 0 |
k3 = 0.2 |
k3 = 0.4 |
k3 = 0.6 |
k3 = 0.8 |
k3 = 1 |
|||
|
SS Beam |
0 |
EBBT (Present) |
2.8034 |
2.8007 |
2.7927 |
2.7794 |
2.7611 |
2.7380 |
2.8462 |
2.8434 |
2.8349 |
2.8209 |
2.8016 |
2.7772 |
|
TBT (Present) |
2.6772 |
2.6740 |
2.6646 |
2.6491 |
2.6278 |
2.6009 |
2.8371 |
2.8343 |
2.8257 |
2.8115 |
2.7920 |
2.7673 |
||
|
TBT [24] |
2.6767 |
2.6748 |
2.6669 |
2.6533 |
2.6337 |
2.6103 |
2.8369 |
2.8349 |
2.8251 |
2.8115 |
2.7919 |
2.7685 |
||
|
RBT |
2.6773 |
2.6746 |
2.6665 |
2.6532 |
2.6347 |
2.6114 |
2.8371 |
2.8343 |
2.8258 |
2.8118 |
2.7925 |
2.7681 |
||
|
0.4 |
EBBT (Present) |
2.7987 |
2.7961 |
2.7882 |
2.7752 |
2.7572 |
2.7345 |
2.8417 |
2.8389 |
2.8304 |
2.8165 |
2.7973 |
2.7729 |
|
|
TBT (Present) |
2.6720 |
2.6687 |
2.6590 |
2.6429 |
2.6207 |
2.5928 |
2.8326 |
2.8297 |
2.8212 |
2.8070 |
2.7874 |
2.7627 |
||
|
TBT [24] |
2.6728 |
2.6689 |
2.6611 |
2.6474 |
2.6279 |
2.6044 |
2.8330 |
2.8291 |
2.8212 |
2.8076 |
2.7880 |
2.7626 |
||
|
RBT |
2.6722 |
2.6694 |
2.6613 |
2.6479 |
2.6293 |
2.6059 |
2.8326 |
2.8298 |
2.8213 |
2.8073 |
2.7880 |
2.7636 |
||
|
0.8 |
EBBT (Present) |
2.7848 |
2.7822 |
2.7745 |
2.7619 |
2.7443 |
2.7221 |
2.8282 |
2.8254 |
2.8171 |
2.8032 |
2.7841 |
2.7599 |
|
|
TBT (Present) |
2.6566 |
2.6532 |
2.6430 |
2.6262 |
2.6032 |
2.5742 |
2.8191 |
2.8162 |
2.8076 |
2.7935 |
2.7739 |
2.7492 |
||
|
TBT [24] |
2.6572 |
2.6533 |
2.6455 |
2.6318 |
2.6123 |
2.5888 |
2.8193 |
2.8154 |
2.8076 |
2.7939 |
2.7744 |
2.7509 |
||
|
RBT |
2.6567 |
2.6540 |
2.6458 |
2.6324 |
2.6138 |
2.5903 |
2.8191 |
2.8162 |
2.8078 |
2.7939 |
2.7747 |
2.7503 |
||
|
1 |
EBBT (Present) |
2.7744 |
2.7718 |
2.7643 |
2.7518 |
2.7345 |
2.7126 |
2.8182 |
2.8154 |
2.8071 |
2.7933 |
2.7742 |
2.7501 |
|
|
|
|
|
|
|
|
|
2.8089 |
2.8061 |
2.7975 |
2.7834 |
2.7639 |
2.7392 |
||
|
TBT (Present) |
2.6450 |
2.6415 |
2.6312 |
2.6141 |
2.5906 |
2.5610 |
2.8095 |
2.8056 |
2.7978 |
2.7841 |
2.7646 |
2.7412 |
||
|
TBT [24] |
2.6455 |
2.6416 |
2.6337 |
2.6201 |
2.6005 |
2.5771 |
2.8089 |
2.8061 |
2.7977 |
2.7838 |
2.7647 |
2.7404 |
||
|
RBT |
2.6452 |
2.6424 |
2.6343 |
2.6208 |
2.6022 |
2.5788 |
|
|
|
|
|
|
||
Note: EBBT: Euler-Bernoulli beam theory; TBT: Timoshenko beam theory; RBT: Refined Beam Theory
Table 3. Dimensionless frequency of different boundary conditions bi-directional functionally graded beam in free vibration
|
Gradient Index, k1 |
Boundary Conditions |
L/h = 5 |
L/h = 20 |
||||||||||
|
Gradient Index, k3 |
Gradient Index, k3 |
||||||||||||
|
k3 = 0 |
k3 = 0.2 |
k3= 0.4 |
k3 = 0.6 |
k3 = 0.8 |
k3 = 1 |
k3 = 0 |
k3 = 0.2 |
k3 = 0.4 |
k3 = 0.6 |
k3 = 0.8 |
k3 = 1 |
||
|
0 |
SS Beam |
2.6773 |
2.6746 |
2.6665 |
2.6532 |
2.6347 |
2.6114 |
2.8371 |
2.8343 |
2.8258 |
2.8118 |
2.7925 |
2.7681 |
|
CC Beam |
5.2361 |
5.2321 |
5.2201 |
5.2003 |
5.1728 |
5.1380 |
6.3512 |
6.3450 |
6.3265 |
6.2958 |
6.2535 |
6.2000 |
|
|
CS Beam |
4.4073 |
4.4029 |
4.3900 |
4.3685 |
4.3389 |
4.3015 |
3.8923 |
3.8890 |
3.8793 |
3.8632 |
3.8409 |
3.8127 |
|
|
CF Beam |
0.9848 |
0.9839 |
0.9811 |
0.9764 |
0.9700 |
0.9619 |
1.0130 |
1.0120 |
1.0090 |
1.0040 |
0.9971 |
0.9884 |
|
|
0.4 |
SS Beam |
2.6722 |
2.6694 |
2.6613 |
2.6479 |
2.6293 |
2.6059 |
2.8326 |
2.8298 |
2.8213 |
2.8073 |
2.7880 |
2.7636 |
|
CC Beam |
5.2403 |
5.2363 |
5.2243 |
5.2045 |
5.1770 |
5.1421 |
6.3573 |
6.3511 |
6.3326 |
6.3019 |
6.2595 |
6.2060 |
|
|
CS Beam |
3.7831 |
3.7800 |
3.7705 |
3.7549 |
3.7333 |
3.7059 |
4.2872 |
4.2830 |
4.2704 |
4.2495 |
4.2207 |
4.1842 |
|
|
CF Beam |
0.8720 |
0.8712 |
0.8687 |
0.8645 |
0.8588 |
0.8516 |
0.8950 |
0.8941 |
0.8915 |
0.8871 |
0.8810 |
0.8733 |
|
|
0.8 |
SS Beam |
2.6567 |
2.6540 |
2.6458 |
2.6324 |
2.6138 |
2.5903 |
2.8191 |
2.8162 |
2.8078 |
2.7939 |
2.7747 |
2.7503 |
|
CC Beam |
5.2531 |
5.2491 |
5.2371 |
5.2172 |
5.1897 |
5.1547 |
6.3758 |
6.3695 |
6.3509 |
6.3202 |
6.2777 |
6.2240 |
|
|
CS Beam |
3.6728 |
3.6697 |
3.6606 |
3.6455 |
3.6245 |
3.5980 |
4.1683 |
4.1642 |
4.1519 |
4.1316 |
4.1036 |
4.0682 |
|
|
CF Beam |
0.7698 |
0.7691 |
0.7669 |
0.7632 |
0.7581 |
0.7518 |
0.7885 |
0.7877 |
0.7854 |
0.7815 |
0.7762 |
0.7694 |
|
|
1 |
SS Beam |
2.6452 |
2.6424 |
2.6343 |
2.6208 |
2.6022 |
2.5788 |
2.8089 |
2.8061 |
2.7977 |
2.7838 |
2.7647 |
2.7404 |
|
CC Beam |
5.2628 |
5.2587 |
5.2467 |
5.2268 |
5.1992 |
5.1642 |
6.3897 |
6.3834 |
6.3648 |
6.3340 |
6.2914 |
6.2376 |
|
|
CS Beam |
3.6169 |
3.6139 |
3.6049 |
3.5901 |
3.5695 |
3.5435 |
4.1090 |
4.1050 |
4.0929 |
K3 = 0.6 |
4.0453 |
4.0103 |
|
|
CF Beam |
0.7225 |
0.7218 |
0.7197 |
0.7163 |
0.7115 |
0.7055 |
0.7394 |
0.7386 |
0.7364 |
2.8118 |
0.7278 |
0.7214 |
|
Note: CC: clamped–clamped; CS: clamped–simply supported, SS: simply supported; CF: clamped-free
In Table 3, the non-dimensional fundamental frequencies (ω) of BDFG beams are presented for various boundary conditions: SS, CC, CS, and CF. The analysis considers two aspect ratios (L/h = 5,20) and several material gradient indices (k1, k3) using the TBT.
The relative difference between the TBT and EBBT frequency results grows with the mode number. For instance, the relative frequency difference for an SS BDFG beam with an aspect ratio of L/h = 5 and material gradient indices of k1 = 0.4 and k3 = 1 is 0.0518%. This relative difference was calculated using the following formula:
Difference $\%=\left[\frac{\bar{\omega}_{E B B T}-\overline{\omega l}_{T B T}}{\overline{\omega l}_{E B B T}}\right] \times 100$ (35)
Table 3 displays the dimensionless frequencies of BDFG beams for various gradient indices (k1, k3), boundary conditions, and slenderness ratios (L/h = 5, 20). The findings reveal that for all examined models, the dimensionless fundamental frequency (ω) drops as the gradient indices (k1, k3) grow. The BDFG beams get more flexible with higher gradient indices because the effective modulus of elasticity goes down, which makes the BDFG beams bend more.
From Figure 3, it is clear how the dimensionless fundamental frequency of a BDFG beam with an aspect ratio L/h = 20 under TBT changes for four different boundary conditions: SS, CC, CS, and CF. The outcomes are shown versus the through-thickness gradient index (k3) and the longitudinal gradient index (k1). The frequency rises monotonically with either gradient index in all circumstances; the highest values are seen when k1 and k3 are both high. The SS condition exhibits the lowest frequencies among the backed scenarios; k3 is more sensitive to a steady increase. The CC beam shows the most frequencies overall because the clamped ends limit bending, and the effect of grading is more noticeable than in SS. Between SS and CC, the CS condition occurs when the reaction is very sensitive to k3 and shows fast gains at low index values and tapers at larger values. The CF beam, which stands for the cantilever case, has the lowest overall frequency values because it has a free end. The frequency surface is almost straight, which means that both gradient indices play a part that is almost the same.
Figure 3. The dimensionless frequency of various aspect ratios (L/h = 20) of the bi-directional functionally graded beam at four boundary conditions in Timoshenko beam theory
Using higher-order beam theory, Table 4 lists the dimensionless frequencies of free vibration beams with various boundary conditions (SS, CS, CC, and CF). The findings demonstrate that the frequencies quickly converge and approach almost constant values after around 6-8 terms as the number of terms in the series solution grows from 2 to 14. This confirms the solution's accuracy and stability. Among the limit conditions, the CC case has the greatest frequencies (about 6.3), whereas the CF case has the lowest values (about 0.73); the (SS ≈ 2.8) and supported (CS ≈ 4.2) cases fall in between. The table also shows that the frequency somewhat drops as the parameter k1 increases, whereas a higher stiffness ratio (ks) somewhat lowers the frequency as well.
Table 4. Dimensionless frequency of different boundaries bi-directional functionally graded beams in free vibration at higher order beam theory
|
Gradient Index, k1 |
No. of Terms (RBT) |
SS Boundary Conditions |
CS Boundary Conditions |
CC Boundary Conditions |
CF Boundary Conditions |
||||
|
k3 = 0.4 |
k3 = 0.8 |
k3 = 0.4 |
k3 = 0.8 |
k3 = 0.4 |
k3 = 0.8 |
k3 = 0.4 |
k3 = 0.8 |
||
|
0 |
2terms |
3.1342 |
3.0972 |
4.0411 |
4.0280 |
6.3617 |
6.2881 |
1.0140 |
1.0021 |
|
4terms |
2.8267 |
2.7933 |
3.9064 |
3.8675 |
6.3333 |
6.2601 |
1.0090 |
0.9972 |
|
|
6terms |
2.8258 |
2.7925 |
3.8896 |
3.8509 |
6.3298 |
6.2567 |
1.0090 |
0.9971 |
|
|
8terms |
2.8258 |
2.7925 |
3.8831 |
3.8446 |
6.3282 |
6.2552 |
1.0090 |
0.9971 |
|
|
10terms |
2.8258 |
2.7925 |
3.8805 |
3.8421 |
6.3273 |
6.2543 |
1.0090 |
0.9971 |
|
|
12terms |
2.8258 |
2.7925 |
3.8796 |
3.8412 |
6.3268 |
6.2538 |
1.0090 |
0.9971 |
|
|
14terms |
2.8258 |
2.7925 |
3.8793 |
3.8409 |
6.3265 |
6.2535 |
1.0090 |
0.9971 |
|
|
0.4 |
2terms |
3.1241 |
3.0871 |
4.3038 |
4.2948 |
6.3667 |
6.2940 |
0.9013 |
0.8907 |
|
4terms |
2.8220 |
2.7887 |
4.2737 |
4.2240 |
6.3394 |
6.2662 |
0.8915 |
0.8810 |
|
|
6terms |
2.8213 |
2.7880 |
4.2718 |
4.2221 |
6.3359 |
6.2628 |
0.8915 |
0.8810 |
|
|
8terms |
2.8213 |
2.7880 |
4.2711 |
4.2214 |
6.3343 |
6.2612 |
0.8915 |
0.8810 |
|
|
10terms |
2.8213 |
2.7880 |
4.2707 |
4.2210 |
6.3334 |
6.2603 |
0.8915 |
0.8810 |
|
|
12terms |
2.8213 |
2.7880 |
4.2705 |
4.2208 |
6.3329 |
6.2598 |
0.8915 |
0.8810 |
|
|
14terms |
2.8213 |
2.7880 |
4.2704 |
4.2207 |
6.3326 |
6.2595 |
0.8915 |
0.8810 |
|
|
0.8 |
2terms |
3.0946 |
3.0579 |
4.1974 |
4.1925 |
6.3822 |
6.3120 |
0.8016 |
0.7922 |
|
4terms |
2.8083 |
2.7751 |
4.1556 |
4.1073 |
6.3580 |
6.2846 |
0.7855 |
0.7762 |
|
|
6terms |
2.8078 |
2.7747 |
4.1535 |
4.1052 |
6.3543 |
6.2810 |
0.7854 |
0.7762 |
|
|
8terms |
2.8078 |
2.7747 |
4.1527 |
4.1044 |
6.3527 |
6.2794 |
0.7854 |
0.7762 |
|
|
10terms |
2.8078 |
2.7747 |
4.1523 |
4.1040 |
6.3518 |
6.2785 |
0.7854 |
0.7762 |
|
|
12terms |
2.8078 |
2.7747 |
4.1520 |
4.1037 |
6.3512 |
6.2780 |
0.7854 |
0.7762 |
|
|
14terms |
2.8078 |
2.7747 |
4.1519 |
4.1036 |
6.3509 |
6.2777 |
0.7854 |
0.7762 |
|
|
1 |
2terms |
3.0730 |
3.0365 |
4.1499 |
4.1466 |
6.3944 |
6.3261 |
0.7525 |
0.7409 |
|
4terms |
2.7980 |
2.7650 |
4.0967 |
4.0491 |
6.3720 |
6.2984 |
0.7329 |
0.7215 |
|
|
6terms |
2.7977 |
2.7647 |
4.0946 |
4.0469 |
6.3683 |
6.2948 |
0.7328 |
0.7214 |
|
|
8terms |
2.7977 |
2.7647 |
4.0937 |
4.0461 |
6.3666 |
6.2931 |
0.7328 |
0.7214 |
|
|
10terms |
2.7977 |
2.7647 |
4.0932 |
4.0456 |
6.3656 |
6.2922 |
0.7328 |
0.7214 |
|
|
12terms |
2.7977 |
2.7647 |
4.0930 |
4.0454 |
6.3651 |
6.2917 |
0.7328 |
0.7214 |
|
|
14terms |
2.7977 |
2.7647 |
4.0929 |
4.0453 |
6.3648 |
6.2914 |
0.7328 |
0.7214 |
|
3.2 Results of forced vibration
Figure 4 shows the fluctuation of the maximum normalised dynamic deflections with the moving load velocity for several beam configurations. Figure 4(a) shows the example of a simply supported homogeneous beam (SS) whose deflection rises with load velocity until it reaches a resonance peak, then declines as the velocity continues to rise. Smaller separations between moving loads (d = L/16) generate greater peak deflections as a result of increased vibration amplification; larger spacing (d = L/2) lowers the response. By contrast, Figure 4(b) depicts the CC-BDFG, where the deflections are notably less than in the SS case due to the greater stiffness given by the clamped borders and the functional grading. Even though resonance peaks are still there, they are very small, and the effect of load spacing is not as clear. The simply supported SS-BDFG in Figure 4(c) has a dynamic response lower than that of the homogeneous SS beam but higher than that of the CC case. Increasing the material gradation's stiffness reduces the peak deflections found in the study by Ying et al. [25] and Rahmani and Pedram [26]. The gradient index is quite important here as well. Figure 4(d) shows the CF-BDFG, which has intermediate deflections between the SS and CC configurations. While resonance effects are clear, they are less than those of the homogeneous SS case. The spacing between moving loads still influences the maximum amplitudes; closer loads boost the vibrations. In Figure 4, results are presented for beams with a length-to-height ratio L/h = 20, stiffness parameters stiffness k1 = 0.5, k3 = 0.5, and rotational speed $\Omega$ = 0. The deflections are evaluated for different spacing distances between the moving loads in the x and z directions.
(a)
(b)
(c)
(d)
Figure 4. Variation of the maximum normalized dynamic deflections of bi-directional functionally graded beams under a moving load
Table 5 provides the maximum normalized mid-span deflections and critical velocities of a BDFG beam for several gradient indices, load spacings, and boundary constraints. The SS's deflections are the greatest of all the instances, reaching values over 2.1 when the load spacing is cut down to L/16; the matching critical velocity drops to 139 m/s. Increasing the load spacing to L/2 lowers the deflection to about 1.40 while maintaining a higher critical velocity of 205 m/s, therefore emphasizing the significant impact of load closeness on resonance. Conversely, the CC has the lowest deflections; its maximum values fall below 0.5 throughout all load spacings, and its critical velocities are much greater, reaching up to 435 m/s at L/2 spacing. This shows how well functional grading and clamped borders can reduce vibration amplitudes. With deflections ranging from 0.89 to 0.92 and critical velocities between 219 and 225 m/s, the CS exhibits intermediate behavior; this suggests that the mixed boundary condition offers higher stiffness than SS but is less effective than CC. The CF records deflections between 7.45 and 8.59, which are the highest among the graded beams, and the corresponding critical speeds are quite low (about 9–45 m/s), indicating the freedom added by the free end.
Table 5. Maximum normal deviations of the functionally graded beam center and corresponding velocities in the binary direction
|
Boundary Conditions |
Distance Between Two Loads |
Power-Law Exponent, k2 |
Power-Law Exponent, k3 |
|
|
0.5 |
||||
|
Max. (w(x,t)/D) |
V$\sigma$ (m/s) |
|||
|
SS Beam |
L/2 |
0.5 |
1.40495040201722 |
205 |
|
L/4 |
1.88514706987099 |
157 |
||
|
L/8 |
2.053737252157936 |
141 |
||
|
L/16 |
2.10240005508007 |
139 |
||
|
CC Beam |
L/2 |
0.30008082302545 |
435 |
|
|
L/4 |
0.423827182234671 |
325 |
||
|
L/8 |
0.479369883007252 |
267 |
||
|
L/16 |
0.495642002703985 |
261 |
||
|
CS Beam |
L/2 |
0.605684967964553 |
297 |
|
|
L/4 |
0.814106298634555 |
229 |
||
|
L/8 |
0.895016938094859 |
219 |
||
|
L/16 |
0.916562309404595 |
225 |
||
|
CF Beam |
L/2 |
7.4558974280466 |
9 |
|
|
L/4 |
7.56710032821635 |
43 |
||
|
L/8 |
8.46473807748718 |
45 |
||
|
L/16 |
8.59844692004846 |
43 |
||
Table 6 shows that the maximum mid-span deflection in the SS functionally graded beam was found to be 1.8461 mm at a load velocity of 80 m/s with a distance between the two loads of d = L/16. The greatest displacement for the CC functionally graded beam under identical circumstances was 0.3107 mm. For the CS, functionally graded beam, too, the greatest mid-span displacement was 0.6322 mm at a moving load velocity of 80 m/s and load spacing of d = L/16. These findings match those found in the 2008 study by Kapuria et al. [24].
Table 6. Maximum mid-span displacement in four types of functionally graded beam at bi-direction
|
Gradient Index, k1 = 0.5 |
Gradient Index, k3 = 0.5 |
||||
|
Beam Theory |
Velocities |
D = L/2 |
D = L/4 |
D = L/8 |
D = L/16 |
|
Maximum Midspan Displacement (mm) |
|||||
|
SS |
Vp= 10 m/s |
0.8836 |
1.1618 |
1.1966 |
1.2409 |
|
Vp = 20 m/s |
0.9366 |
1.1236 |
1.2429 |
1.2652 |
|
|
Vp = 40 m/s |
0.8888 |
1.2791 |
1.2626 |
1.3636 |
|
|
Vp = 80 m/s |
0.9362 |
1.4710 |
1.7650 |
1.8461 |
|
|
CC |
Vp = 10 m/s |
0.1541 |
0.2566 |
0.2906 |
0.2996 |
|
Vp = 20 m/s |
0.1542 |
0.2579 |
0.2905 |
0.3012 |
|
|
Vp = 40 m/s |
0.1565 |
0.2598 |
0.2953 |
0.3027 |
|
|
Vp = 80 m/s |
0.1653 |
0.2822 |
0.2997 |
0.3107 |
|
|
CS |
Vp = 10 m/s |
0.3429 |
0.5034 |
0.5500 |
0.5637 |
|
Vp = 20 m/s |
0.3447 |
0.5078 |
0.5532 |
0.5643 |
|
|
Vp = 40 m/s |
0.3480 |
0.5272 |
0.5618 |
0.5655 |
|
|
Vp = 80 m/s |
0.4185 |
0.5424 |
0.6071 |
0.6322 |
|
|
CF |
Vp = 10 m/s |
6.4066 |
6.3277 |
6.8863 |
7.1192 |
|
Vp = 20 m/s |
6.0755 |
6.4785 |
7.0488 |
7.3958 |
|
|
Vp = 40 m/s |
5.7947 |
7.5123 |
8.3303 |
8.5726 |
|
|
Vp = 80 m/s |
5.4218 |
5.8196 |
5.5231 |
4.9776 |
|
Two sequential moving harmonic loads at varied spacing are applied to an SS BDFG beam with k1 = 0.5, k3 = 0.5, L/h = 20, and Ω = 0 in Figure 5, which shows the mid-span displacement time histories. Figures 5(a)-(d) match load distances of d = L/2, d = L/4, d = L/8, and d = L/16, respectively. The findings reveal that displacement responses rely heavily on load spacing as well as load velocity. The two loads behave practically independently at larger spacing (d = L/2), generating smoother displacement curves. Interaction effects get more apparent as the spacing shrinks (d = L/16), which results in bigger peak displacements and more vigorous oscillations, especially at faster load velocities (Vp = 80 m/s). This suggests that the beam's dynamic response is increased by tighter load spacing, which agrees with the research by Kapuria et al. [24].
(a)
(b)
(c)
(d)
Figure 5. Time history of the mid-span displacements of a simply supported bi-directional functionally graded beam for k1 = 0.5, k3 = 0.5, L/h = 20, Ω = 0, a) d = L/2, b) d = L/4, c) d = L/8, d) d = L/16
For various impact positions, Figure 6 displays the time history of mid-span displacements of a CC BDFG beam under various projectile speeds. The beam shows the greatest displacement response with obvious oscillations in Figure 6(a), where the impact happens at the mid-span (d = L/2). This is because this is where the beam is the most flexible. The displacement amplitude drops as the point of impact moves closer to the supports. Though greater velocities still generate noticeable peaks, the mid-span displacement in Figure 6(b) (d = L/4) is smaller, and the oscillations are smoother than in the mid-span impact. The displacement continues to drop in Figure 6(c) (d = L/8), and the response turns more sinusoidal, dominated by the fundamental bending mode. At last, in Figure 6(d) (d = L/16), where the influence is quite near the support, the mid-span displacement is the lowest, and the vibration response is smooth with very little higher-mode effects owing to the strong restraining action of the supports. Figure 7 shows the mid-span displacement time histories for a simply supported functionally graded beam (L/h = 20, d = L/2, P0 = 100 kN) under moving loads with Vp = 10, 20, and 40 m/s with excitation Ω = 0 and 40 rad/s. In Figure 6(a) (k1 = 1, k3 = 0; axial gradient, the responses at Ω = 0 formed smooth, single-hump traces whose peak amplitudes increased with load velocity; at Ω = 40, the response became strongly oscillatory, and its amplitude again grew with Vp. In Figure 6(b) (k1 = 0, k3 = 1; thickness gradient), the same patterns were seen: higher Vp made bigger peaks at Ω = 0 and bigger wavy envelopes at Ω = 40. The magnitudes were almost the same as in Figure 6(a), which shows that the stiffness of the material is about the same for axial and thickness gradients. In Figure 6(c) (k1 = 1 k3 = 1; bi-directional gradient), the displacement histories followed the same qualitative behavior but with slightly smaller maxima than in Figures 6(a) and (b), both for the quasi-static-like case (Ω = 0) and for the harmonic case (Ω = 40). Increasing the moving-load velocity generally enhanced the mid-span response. The harmonic excitation at Ω = 40 showed clear periodic oscillations, axial and thickness gradations gave almost identical amplitudes, but the bi-directional gradation produced somewhat diminished deflections, which supports the findings of Alshorbagy et al. [27].
(a)
(b)
(c)
(d)
Figure 6. Time history of the mid-span1displacements of a clamped–clamped bi-directional functionally graded beam for k1 = 0.5, k3 = 0.5, L/ h = 20, Ω = 0; a) d = L/2, b) d = L/4, c) d = L/8, and d) d = L/16
Figure 8 shows the normalised dynamic deflection of a simply supported functionally graded beam under a moving load speed of Vp = 20 m/s, with a span-to-thickness ratio of h/L = 20, varying distances between two loads 1 (d = L/2, L/4, L/8, L/16), and several power-law indices (k2 = 0, 2, 4, 6, 8, 10) at zero excitation frequency (Ω = 0). The beam showed the greatest deflection in Figure 8(a) when the load separation was d = L/2. As the power-law index rose, the deflection progressively dropped, suggesting that more ceramic material lowered flexibility. Corresponding to d = L/4, Figure 8(b) exhibited lesser deflections than d = L/2, yet the same declining pattern with growing k2. With d = L/8 in Figure 8(c), the deflections were lower, therefore reflecting the lower impact of moving loads set closer together. The dwell time of the moving harmonic load on the functionally graded beam seems to grow as its velocity value falls, which results in more induced vibration cycles. Şimşek and Cansiz [28] agree with this. A reduction in load speed produces less dynamic response than an increase does. In Figure 8(d), with d = L/16, the least deflections were shown. The influence of the power-law index was the clearest as larger k2 values caused a noticeable decrease in response, confirming the findings of Stojanović et al. [30].
(a)
(b)
(c)
(a)
(b)
(c)
(d)
Figure 8. The normalized dynamic deflection of a simply supported functionally graded beam when various velocities of moving loads (Vp = 20, m/sec) and distance between two loads, h/L = 20, with different power law index (k2 = 0, 2, 4, 6, 8, 10) and excitation frequency Ω = 0
Using Euler–Bernoulli, Timoshenko, and higher-order shear deformation theories, this study provided a thorough examination of the free and forced vibration response of BDFG beams under repeated moving harmonic loading. Using Lagrange's equations and the Newmark-β approach, the dynamic response under different boundary conditions, gradient indices, load velocities, excitation frequencies, and load spacings was captured.
Increasing gradient indices lowers the basic frequency, hence enhancing the beam's flexibility. The critical load speed is greatly affected by boundary conditions and the spacing between moving loads; thicker beams show a significant impact from shear deformation. Furthermore, the stress and displacement distributions of BDFG beams vary noticeably from those of homogeneous and one-directional functionally graded beams, therefore highlighting the need for bi-directional gradation for customizing structural performance.
This study's main contributions are:
(i) Comparing the three beam theories' ability to capture BDFG beams' forced and free vibrations.
(ii) A thorough parametric study examining how dynamic loading characteristics, boundary conditions, geometric ratios, and gradient indices affect performance.
(iii) New understanding of how well BDFG beams perform when they are subjected to moving harmonic loads one after another is an area that hasn't been studied much before.
Notwithstanding these improvements, some restrictions still exist. The research uses a linear theoretical framework and ignores material damping, which could affect dynamic responses in real-world situations. Furthermore, since no experimental confirmation was conducted, the practical relevance of the findings to actual structures is restricted.
Future studies should overcome these limits by broadening the analysis to include nonlinear vibration regimes, adding damping models, and confirming predictions via empirical testing. Research on thermal loading, temperature-dependent material properties, and three-dimensional modeling will help us to better understand BDFG structures and extend their use in modern engineering systems.
[1] Mohammadi, M., Rajabi, M., Ghadiri, M. (2021). Functionally graded materials (FGMs): A review of classifications, fabrication methods and their applications. Processing and Application of Ceramics, 15(4): 319-343. https://doi.org/10.2298/PAC2104319M
[2] Ninh, V.T.A. (2023). Dynamics of functionally graded beams carrying a moving load with influence of porosities and partial foundation support. Vietnam Journal of Science and Technology, 61(4): 692-707. https://doi.org/10.15625/2525-2518/17507
[3] Gottron, J., Harries, K.A., Xu, Q. (2014). Creep behaviour of bamboo. Construction and Building Materials, 66: 79-88. https://doi.org/10.1016/j.conbuildmat.2014.05.024
[4] Nguyen, D.K., Vu, A.N. T., Le, N.A.T., Pham, V.N. (2020). Dynamic behavior of a bidirectional functionally graded sandwich beam under nonuniform motion of a moving load. Shock and Vibration. https://doi.org/10.1155/2020/8854076
[5] Bohidar, S.K., Sharma, R., Mishra, P.R. (2014). Functionally graded materials: A critical review. International Journal of Research, 1(7): 289-301.
[6] Şimşek, M., Al-Shujairi, M. (2017). Static, free and forced vibration of functionally graded (FG) sandwich beams excited by two successive moving harmonic loads. Composites Part B: Engineering, 108: 18-34. https://doi.org/10.1016/j.compositesb.2016.09.098
[7] Şimşek, M. (2015). Bi-directional functionally graded materials (BDFGMs) for free and forced vibration of Timoshenko beams with various boundary conditions. Composite Structures, 133: 968-978. https://doi.org/10.1016/j.compstruct.2015.08.021
[8] Deng, H., Cheng, W. (2016). Dynamic characteristics analysis of bi-directional functionally graded Timoshenko beams. Composite Structures, 141: 253-263. https://doi.org/10.1016/j.compstruct.2016.01.051
[9] Nejad, M.Z., Hadi, A., Rastgoo, A. (2016). Buckling analysis of arbitrary two-directional functionally graded Euler–Bernoulli nano-beams based on nonlocal elasticity theory. International Journal of Engineering Science, 103: 1-10. https://doi.org/10.1016/j.ijengsci.2016.03.001
[10] Nejad, M.Z., Hadi, A. (2016). Non-local analysis of free vibration of bi-directional functionally graded Euler-Bernoulli nano-beams. International Journal of Engineering Science, 105: 1-11. https://doi.org/10.1016/j.ijengsci.2016.04.011
[11] Nejad, M.Z., Hadi, A. (2016). Eringen’s non-local elasticity theory for bending analysis of bi-directional functionally graded Euler–Bernoulli nano-beams. International Journal of Engineering Science, 106: 1-9. https://doi.org/10.1016/j.ijengsci.2016.05.005
[12] Wang, Z., Wang, X., Xu, G., Cheng, S., Zeng, T. (2016). Free vibration of two-directional functionally graded beams. Composite Structures, 135: 191-198. https://doi.org/10.1016/j.compstruct.2015.09.013
[13] Shafiei, N., Kazemi, M. (2017). Buckling analysis on the bi-dimensional functionally graded porous tapered nano-/micro-scale beams. Aerospace Science and Technology, 66: 1-11. https://doi.org/10.1016/j.ast.2017.02.019
[14] Shafiei, N., Mirjavadi, S.S., Mohasel Afshari, B., Rabby, S., Kazemi, M. (2017). Vibration of two-dimensional imperfect functionally graded (2D-FG) porous nano-/micro-beams. Computational Methods in Applied Mechanics and Engineering, 322: 615-632. https://doi.org/10.1016/j.cma.2017.05.007
[15] Van Do, T., Nguyen, D.K., Duc, N.D., Doan, D.H., Bui, T.Q. (2017). Analysis of bi-directional functionally graded plates by FEM and a new third-order shear deformation plate theory. Thin-Walled Structures, 119: 687-699. https://doi.org/10.1016/j.tws.2017.07.022
[16] Nguyen, D.K., Nguyen, Q.H., Tran, T.T., Bui, V.T. (2017). Vibration of bi-dimensional functionally graded Timoshenko beams excited by a moving load. Acta Mechanica, 228(1): 141-155. https://doi.org/10.1007/s00707-016-1705-3
[17] Trinh, L.C., Vo, T.P., Thai, H.T., Nguyen, T.K. (2018). Size-dependent vibration of bi-directional functionally graded microbeams with arbitrary boundary conditions. Composite Structures Part B, 134: 225-245. https://doi.org/10.1016/j.compositesb.2017.09.054
[18] Yang, T., Tang, Y., Li, Q., Yang, X.D. (2018). Nonlinear bending, buckling and vibration of bi-directional functionally graded nanobeams. Composite Structures, 204: 313-319. https://doi.org/10.1016/j.compstruct.2018.07.045
[19] Taima, M.S., Shehab, M.B., El-Sayed, T.A., Friswell, M.I. (2023). Comparative study on free vibration analysis of rotating bi-directional functionally graded beams using multiple beam theories with uncertainty considerations. Scientific Reports, 13(1): 17917. https://doi.org/10.1038/s41598-023-44411-0
[20] Huang, X., Zhang, L., Ge, R. (2023). Investigation of the free vibrations of bi-directional functionally graded material beams using differential quadrature method. Journal of Vibration and Control, 30(21-22): 4828-4851. https://doi.org/10.1177/10775463231214314
[21] Belabed, Z., Tounsi, A., Bousahla, A.A., Tounsi, A., Bourada, M., Al-Osta, M.A. (2024). Free vibration analysis of bi-directional functionally graded beams using a simple and efficient finite element model. Structural Engineering and Mechanics, 90(3): 233-252. https://doi.org/10.12989/sem.2024.90.3.233
[22] Yang, J., Chen, Y. (2008). Free vibration and buckling analyses of functionally graded beams with edge cracks. Composite Structures, 83(1): 48-60. https://doi.org/10.1016/j.compstruct.2007.03.006
[23] Lohar, H., Mitra, A. (2025). Free vibration analysis of pre-deformed bi-directionally graded beam supported on nonlinear elastic foundation. Journal of the Institution of Engineers (India): Series C, 106(1): 157-169. https://doi.org/10.1007/s40032-024-01134-z
[24] Kapuria, S., Bhattacharyya, M., Kumar, A.N. (2008). Bending and free vibration response of layered functionally graded beams: A theoretical model and its experimental validation. Composite Structures, 82(3): 390-402. https://doi.org/10.1016/j.compstruct.2007.01.019
[25] Ying, J., Lü, C.F., Chen, W.Q. (2008). Two-dimensional elasticity solutions for functionally graded beams resting on elastic foundations. Composite Structures, 84(3): 209-219. https://doi.org/10.1016/j.compstruct.2007.07.004
[26] Rahmani, O., Pedram, O. (2014). Analysis and modeling the size effect on vibration of functionally graded nanobeams based on nonlocal Timoshenko beam theory. International Journal of Engineering Science, 77: 55-70. https://doi.org/10.1016/j.ijengsci.2013.12.003
[27] Alshorbagy, A.E., Eltaher, M.A., Mahmoud, F.F. (2011). Free vibration characteristics of a functionally graded beam by finite element method. Applied Mathematical Modelling, 35(1): 412-425. https://doi.org/10.1016/j.apm.2010.07.006
[28] Şimşek, M., Cansiz, S. (2012). Dynamics of elastically connected double-functionally graded beam systems with different boundary conditions under action of a moving harmonic load. Composite Structures, 94(9): 2861-2878. https://doi.org/10.1016/j.compstruct.2012.03.016
[29] Şimşek, M. (2010). Vibration analysis of a functionally graded beam under a moving mass by using different beam theories. Composite Structures, 92(4): 904-917. https://doi.org/10.1016/j.compstruct.2009.09.030
[30] Stojanović, V., Dimitrov, L., Tomov, P., Li, D., Nikolić, V. (2024). Dynamic stability analysis of a coupled moving bogie system. Facta Universitatis, Series: Mechanical Engineering, 22(4): 773-785. https://doi.org/10.22190/FUME241003045S