Finite Element Modeling of Dynamic Deformation of Bidirectional Functionally Graded Beams
- Authors: Vasserman I.N.1, Glot I.O.1, Shardakov I.N.1
-
Affiliations:
- Institute of Continuous Media Mechanics, Ural Branch of the Russian Academy of Sciences (UB RAS)
- Issue: Vol 22, No 4 (2026)
- Pages: 317-330
- Section: Analytical and numerical methods of analysis of structures
- URL: https://journals.rudn.ru/structural-mechanics/article/view/53141
- DOI: https://doi.org/10.22363/1815-5235-2026-22-4-317-330
- EDN: https://elibrary.ru/LKNCQZ
- ID: 53141
Cite item
Abstract
This article presents a finite element algorithm for the numerical solution of boundary value problems for the deformation state of elastic beams made of functionally graded materials. The presented algorithm assumes that the density and elastic properties of the beam depend on the coordinates of the cross-section and on the coordinate along the longitudinal axis. The presented formulation allows to solve static and dynamic problems. The Euler-Bernoulli model and the third-order Reddy model are used to describe the deformation state of the beam. The beam finite element utilizes a Hermitian approximation. Using Hamilton's variational principle, the beam equation of motion is obtained in matrix form. The numerical implementation of the resulting finite element algorithm is carried out in Python using the Numpy and Scipy libraries. As an example, a numerical analysis of the solution to the problem of the natural frequencies of vibrations of a bidirectional functionally graded cantilever beam is performed. The numerical convergence and reliability of the obtained solutions are demonstrated. The possibility of coupling of the longitudinal and transverse vibrations of a beam is demonstrated for certain distributions of the physical properties of the material in the beam’s cross-section. It is shown that distributing the material properties along the longitudinal axis improves the beam’s dynamic properties.
Full Text
1. Introduction Functionally graded materials (FGMs) are a class of composites whose properties vary continuously and smoothly in at least one spatial direction. These property variations may be associated with composition, porosity, or microstructure. Such materials allow to optimally combine component properties and minimize stress concentration at the interface between different materials [1; 2]. Unlike layered composites, they are more resistant to crack propagation and delamination [3]. The concept of functionally graded materials was first proposed as early as the 1970s [4]. However, this idea had not yet been implemented at that time due to the lack of appropriate technologies. It was carried into effect in the following decade in Japan as part of an aerospace project [5; 6], which required the creation of a thermal barrier in a 10 mm thick structure wall. Currently, functionally graded materials are attracting increasing attention from researchers around the world. Examples of applications of functionally graded materials can be found in [7]. These include prosthetics that combine mechanical strength in the inner layers with biocompatibility in the outer layers [8], functionally graded blades for turbomachinery and wind turbines [9], reinforcement of mine tunnel walls [10], piezoelectric elements for structural vibration monitoring [11], microsensors and microactuators in nano- and microelectromechanical systems [12]. Recently, there has been growing interest in the use of functionally graded materials in building structures [13]. In particular, much attention is given to the development of structures made of functionally graded concrete. Rational distribution of cement content in a structure allows to reduce its production cost and energy consumption [14], minimize residual stresses, and prevent cracking during the setting process [15]. Additive manufacturing and 3D printing, which are often used in the creation of graded materials [16; 17], are well-suited for concrete structures [18-20]. The heterogeneity of functionally graded materials used in building structures is often associated with porosity [21]. Structures made of such materials have a good weight-to-strength ratio, as well as good vibration- and noise-absorbing properties. On the other hand, theories of functionally graded porous materials can be used in modeling distributed defects. Beams of short to medium length with rectangular cross-sections are typical elements of structures. Treating them as three-dimensional spatial structures results in unjustified computational costs. The earliest and most widely used theory of beams is the classical Euler - Bernoulli theory, in which the cross-section is assumed to be perpendicular to the neutral axis of the beam. This theory does not account for shear in the cross-section; it is suitable for long, slender beams with a high ratio of the longitudinal dimension to the transverse dimensions. For beams of short and medium length, as well as for cases where the assumed deformation shape of the beam neutral axis is comparable to the dimensions of its cross-section, this theory is insufficient. The first-order beam theory (Timoshenko) allows to analyze such structures. However, it requires knowledge of an additional parameter - the shear correction factor - which cannot be determined in terms of this theory. Furthermore, since the shear strain within the cross-section is assumed to be constant, the boundary conditions for the absence of shear stresses at the boundary of the cross-section are not satisfied [22]. All of the aforementioned problems are particularly relevant for heterogeneous beams, including functionally graded beams. Shear deformation theories of higher-order are free of many shortcomings. In these theories, the distribution of displacements across the cross-section is nonlinear and ensures that shear stresses at the free boundary are zero [23-25]. One such theory is Reddy’s third-order theory [26]. The presented review highlights the growing need to develop and implement effective numerical algorithms for solving static and dynamic boundary value problems involving the deformation of structures made of functionally graded materials. The capabilities of well-known software tools such as ANSYS and COMSOL for solving such problems are quite limited. Most articles on the dynamics of the considered structures focus on unidirectional functionally graded beams. Bidirectional functionally graded beams have been studied in a small number of papers [27-30]. These studies approximate displacements and rotation angles using series of functions that satisfy the kinematic boundary conditions in advance. This approach ensures good accuracy, but significantly reduces the versatility of the numerical method. This is particularly important when the beam is considered not as an isolated object, but as part of a more complex structure. The finite element method provides a single approach for beams with different boundary conditions. This article presents a variational formulation and implements a finite element method for beams made of functionally graded materials. A heterogeneous beam is considered, the physical properties of which can vary both within the cross-section and along the longitudinal axis. It will be referred to as a bidirectional functionally graded beam. The Euler - Bernoulli model and Reddy’s third-order model [26] are used to describe the deformed state of this beam. Hermitian approximation is used in the beam finite element. The numerical implementation of the resulting finite-element algorithm was carried out in Python. A numerical analysis of the solution to the problem of natural modes and frequencies of a bidirectional functionally graded beam was performed. The numerical convergence of the obtained solutions and their reliability were demonstrated. The possibility of coupling of the longitudinal and transverse vibrations of the beam is shown for certain distributions of physical properties of the material across the beam section. 2. Methods Beams of rectangular cross-section are widely used in construction (Figure 1). The figure shows the coordinate axes and the characteristic dimensions of the structure. Figure 1. Bidirectional functionally graded beam of rectangular cross-section S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. Here L is the beam length; h and b are the height and width. Directions of the coordinate axes: x is directed along the longitudinal axis of the beam; z is perpendicular to the beam axis and lies in the plane of bending; y is perpendicular to the pane of bending. 2.1. Properties of Functionally Graded Materials A typical functionally graded material is a composite consisting of particles with a specified continuous distribution of the volume fractions of the constituent phases. In the case of isotropic phases, the material is isotropic, but heterogeneous. Among the models that describe the effective properties of a material based on the properties of its components, the most common and simplest is the Voigt model (the rule of mixtures). Pf P1V1 P2V2 , (1) where Pf is the effective property of the material; Pi is the property of the pure i-th component of the material; Vi is the volume fraction of the pure i-th component of the material. The relationship V1 V2 1 holds. Young’s modulus E , shear modulus G and mass density can serve as property P . There are also a number of other models for describing the effective properties of a material, a review of which can be found in [31]. For moderate porosity, the theory of mixtures can also be used: Pf Pf 0 1 V0 , (2) where Pf 0 is the effective property of the material without pores, calculated according to (1), V0 is the volume fraction of pores. More complex relationships between material properties and porosity, taking into account the pore structure, can be found in [13, Table 3] and in [21]. 2.2. Beam Models In the classical Euler-Bernoulli theory, the distribution of displacements across the cross-section takes the form: ux x,z u0 x zw,0x x ; (3) uz x, z w0 x . (4) Here ux and uz are the longitudinal and transverse displacements of an arbitrary point of the beam; u0 and w0 are the longitudinal and transverse displacements of a point on the central axis of the beam. The deformation energy and kinetic energy for them are expressed, according to [32], as 1L 0 2 0 0 0 2 U 2 0 D11 u,x 2D12 u,x w,xx D22 w,xx dx (5) and K 1 L I11 u&0 2 I11 w& 2 2I12 u&0 w&,0x I22 w&,0x 2 dx , (6) 0 where D11, D12, D22, I11, I12, I22 are the integral parameters of the cross-section, calculated as D11 x E x,z dA; D12 x E x, z zdA; D22 x E x,z z2dA, A A A 2 (7) 2 I11 x x,z dA; I12 x x,z zdA; I22 x x,z z dA. (8) A A A It is evident that, unlike homogeneous beams, when the material is distributed asymmetrically along the height, terms that couple longitudinal and transverse vibrations appear. In higher-order theories, the distribution of displacements across the section can be represented as follows [22-26]: ux x, z u0 x zw,0x x z 0 x , uz x, z w0 x . (10) Here 0 x w,0x x 0 x (11) is the transverse shear strain of a point on the neutral axis of the beam; 0 is the total flexural rotation of the cross-section on the neutral axis; z is the correction function. Angles 0 and 0 are shown in Figure 2. Function z depends on the cross-section and is selected in a way such that the condition of zero shear stresses on the boundary is satisfied. Taking z 0 results in the Euler - Bernoulli theory, z z results in the Timoshenko theory. For Reddy’s third-order theory [26], the correction function has the form: z z 4z2 . (12) 1 3h2 (9) Figure 2. Deformation of the cross-section: 1 - center line; 2 - normal plane; 3 - deformed cross-section; 4 - plane tangent to the deformed cross-section at the point of its intersection with the center line. S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. The deformation energy of a beam with a rectangular cross-section can be expressed as [24]: 1 L UR Axx u,0x 2 2Bxx u,0x w,0xx 2Exx u,0x ,0x Dxx w,0xx 2 2 0 2 2 2Fxx w,0xx ,0x Hxx ,0x Axz 0 dx, accordingly, the kinetic energy is: R 1 L IA u&0 2 2IB u&0 w&,0x 2IE u&0 &0 ID w&,x 2 K 2 0 (13) 2IF w&,0x &0 IH &0 2 IA w& 0 2 dx , (14) where Axx, Bxx, Dxx, IA, IB, ID are the integral parameters of the section, similar to parameters D11, D12, D22, I11, I12, I22 for Bernoulli beams: Axx x E x,z dA; Bxx x E x,z zdA; Dxx x E x, z z2dA; A A A (15) IA x x,z dA; IB x x, z zdA; ID x x,z z2dA. (16) A A A Exx, Fxx, Hxx , IE , IF , IH are the integral parameters of the section, containing correction function z : Exx x E x,z z dA; Fxx x E x,z z zdA; Hxx x E x, z z 2dA; A A A (17) IE x x,z z dA; IF x x,z z zdA; IH x x,z z 2dA, (18) A A A Axz is the integral parameter of the section, associated with the contribution of shear strain, containing the derivative of the correction function with respect to z : Axz x G x, z ,z z 2dA. (19) A In expressions (5, 6, 13, 14), terms D12, I12, Bxx, IB, Exx, IE couple the longitudinal and transverse vibrations. In the case of a material distribution across the section that is symmetric about the central axis, these terms are zero, and the corresponding vibrations are uncoupled. With an asymmetric distribution of material properties, the distribution of axial stresses will also be asymmetric, leading to the coupling of bending and tension in the beam. When considering bending not about the geometric center of the cross-section, but about the physical center zc , for which E z z zc dA 0, then bending and tension A can be decoupled for static problems. An approach that allows decoupling of bending and tension for higher-order bending theory is proposed in [33]. However, for dynamic problems, it is possible to simultaneously reduce D12 and I12, Bxx and IB , Exx and IE to zero in this manner only if the ratio of Young’s modulus to density is constant across the section. In a general case, the interaction between bending and tension will occur regardless of the choice of the centerline. This interaction should be distinguished from the effect of longitudinal load on the natural frequency of transverse vibration, which is closely related to Euler instability. 2.3. Finite Element Approximation and Numerical Solution The finite element method is used to numerically solve the free vibration problem. The element diagram is shown in Figure 3. Each node is described by ui,wi,wi w / x i . In Reddy’s model either the total bending rotation at the node i , or the transverse shear strain of a point on the neutral axis i are added to the nodal unknowns. Let these variants of Reddy’s model be called - variant and - variant, respectively. For the finite-element approximation of longi- tudinal displacements and total bending rotations of the cross-section, the Lagrange approximation is used; for transverse displacements, the Hermitian approxi- mation is used. u0 u1N1 u2N2, (20) Figure 3. Beam element with two nodes S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. 0 1N1 2N2 , or 0 1N1 2N2 , (21) Here w0 w1H1 w2H3 w1 lˆH2 w2 lˆH4 . (22) N1 1 ; N2 1 , (23) H1 2 3 3 ; H2 1 2 3 ; H3 2 3 3 ; H4 1 2 3 . (24) Using the expressions for deformation energy and kinetic energy (4, 5, 11, 12), Hamilton’s principle, and the finite-element approximation, the equation of motion can be written in matrix form: Mq&& t Kq t F t , (25) where K is the stiffness matrix; M is the mass matrix; q t is the vector of nodal displacements, F t is the vector of external forces. When considering free vibrations, the vector of external forces is equal to zero, and the vector of nodal displacements can be represented as q t q~cos t , (26) where is the natural angular frequency of vibration of the beam, is the phase, vector q~ does not depend on time. As a result, the frequency equation can be obtained in the form K 2M q~ 0 . The local stiffness and mass matrices of the element are calculated as follows: (27) ke ΒTDBdL, L (28) me NTI NdL. (29) L Assembly of the global K and M matrices from the local ke and me matrices is performed in the standard way [34]. For the Bernoulli model: D11 D12 D ; (30) D12 D22 I11 0 I12 ; (31) I 0 I11 0 I12 0 I22 N1 0 0 N2 0 0 ˆH2 0 H3 lˆH4 ; (32) N 0 H1 l 0 H1 lˆH2 0 H3 lˆH4 N1 0 0 N2 0 0 B lˆH2 H3 lˆH4 . (33) H1 For higher-order models: Axx Bxx Exx 0 B D xx Dxx Fxx 0 ; (34) Exx Fxx Hxx 0 0 0 0 Axz IA IB IE 0 IB ID IF 0 . (35) I IE IF IH 0 0 0 0 IA Matrices N and B for the - variant are calculated as N1 0 0 0 N2 0 0 0 0 H1 lˆH2 0 0 H3 lˆH4 0 ; (36) N ˆH2 N1 0 H3 lˆH4 N2 0 H1 l ˆH2 0 0 H3 lˆH4 0 0 H1 l N1 0 0 0 N2 0 0 0 0 H1 lˆH2 0 0 H3 lˆH4 0 , (37) B ˆH2 N1 0 H3 lˆH4 N2 0 H1 l ˆH2 N1 0 H3 lˆH4 N2 0 H1 l and for the - variant: N1 0 N 0 0 0 H1 0 H1 0 lˆH2 0 lˆH2 0 0 N1 0 N2 0 0 0 0 H3 0 H3 0 lˆH4 0 lˆH4 0 0 ; N2 0 (38) N1 0 B 0 0 0 H1 0 0 0 lˆH2 0 0 0 0 N1 N1 N2 0 0 0 0 H3 0 0 0 lˆH4 0 0 0 0 . N2 N2 (39) 2 L N H H Here lˆ , N , H , H , where L is the length of the element. 2 x x x2 A program based on the finite element method was implemented in Python (version 2.7.13) using the NumPy (version 1.11.3) and SciPy (version 0.18.11) libraries. Programs from [34] were used as a prototype. The general structure of the programs for assembling finite element matrices was taken from this book, with Matlab matrix operations replaced by the corresponding Numpy functions. The transition from Matlab to Python is due to the free availability of the latter. The scipy.linalg.eig function with default arguments was used as the eigenvalue solver. Integration of the finite element matrices (28, 29) was performed using the Gauss scheme with a number of points np 2. The cross-sectional characteristics (7, 8, 15-19) were calculated using the Gauss - Lagrange integration scheme (an equivalent of the Gauss scheme that includes points at the ends of the interval) with a number of points nh 11. 3. Results and Discussion The examples below show the results of solving the problem of the natural vibrations of a cantilever beam made of a functionally graded material whose properties vary in both transverse and longitudinal directions. The beam has a rectangular cross-section. In general, the spatial distribution of the beam material properties is described by the following equations: ~ ~ x x E x,z fE z exp L ; x, z f z exp L , (40) ~ where is the parameter characterizing the intensity of variation of material properties along the beam axis; fE z E 0,z and f 0, z are the functions characterizing the variation of the Young’s modulus and the mass density in the transverse direction; L is the beam length. Here, the exponential decrease in the Young’s modulus and mass density in the x-direction from the fixed end to the free end given by (2) can be attributed to an increase in porosity in that direction according to the law ~ x V0 1 exp L , (41) ~ when 0 . In the z-direction, the change in material properties is due to the varying content of components and is asymmetrical relative to the center line of the beam. Functions fE , f were calculated according to (1), where the content of the components varied according to a power law: k z 1 V2 , V1 1 V2, (42) h 2 where k is a non-negative exponent. On the lower surface, the material consists only of the first component; on the upper surface, it consists only of the second. For a beam analyzed using Reddy’s model, Poisson’s ratio was assumed to be constant, and the shear modulus was calculated as E G . 2 1 Dimensionless frequencies were obtained (43) L2, (44) where is the natural angular frequency of vibration of the beam, L is the beam length, and B D, are the integral parameters calculated according to formulas h2 h2 B b f dz , D b z f dz2 E , (45) h2 h2 where b and h are the width and height of the beam, respectively. In the considered examples, material constants from [32] were used to describe the beam material, corresponding to the properties of steel at the lower surface and alumina at the upper surface: El 210 GPa, l 7800kg/m3, Eu 390GPa, u 3980kg/m3. Poisson’s ratio was taken as 0.3, and the exponent k in relation (42) is equal to 3. The results of the calculation of the reduced natural frequencies of the beam obtained using different models are presented below. Table 1 contains the natural frequencies for a beam that is heterogeneous in two directions ( = -1), obtained without taking into account the cross-coupling of longitudinal and transverse vibrations. The numerical solution based on the simplified Bernoulli model (which does not account for the rotational inertia of the cross-section relative to the beam neutral axis) virtually coincides with the result of the analytical solution [30]. Accounting for rotational inertia in the case of short beams leads to a significant decrease in natural frequencies, especially for higher modes. Table 1. Dimensionless natural frequencies of transverse vibrations (Bernoulli model) = -1 Mode Analytical solution [30] Finite element solution using a simplified model Taking into account the rotational inertia L/h=5 Taking into account the rotational inertia L/h=10 Taking into account the rotational inertia L/h=20 1 4.7349 4.7349 4.6924 4.7242 4.7323 2 24.2018 24.2018 23.049 23.8984 24.1250 3 63.8645 63.8645 57.2544 61.9988 63.3825 4 123.0979 123.0981 102.2526 116.6982 121.4002 S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. Asymmetric distribution of the physical properties of the material relative to the centerline leads to coupling of longitudinal and transverse vibrations. Table 2 lists the natural frequencies of the functionally graded beam. Columns 1 and 2 show the values obtained using the model that neglects the interdependence of flexural and longitudinal vibrations, while Column 3 shows the values obtained using the model that accounts for this factor. Column 4 compares the amplitudes of longitudinal and transverse vibrations calculated using the full Bernoulli model. The table shows that accounting for the interdependence of longitudinal and flexural vibrations in an heterogeneous beam leads to a decrease in the natural frequencies, and this factor is more significant for higher modes. Table 2. Interaction of longitudinal and transverse vibrations (Bernoulli model) = -1 Mode Transverse vibrations taking into account rotational inertia Longitudinal vibrations Full Bernoulli model |w|max 1 4.6924 4.6143 0.016 2 23.049 22.5405 0.09 3 31.8648 32.0069 4.6 4 57.2544 55.6523 0.12 5 81.0948 81.7920 4.4 6 102.2526 98.5168 0.24 S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. The analysis of the natural vibrations of a heterogeneous beam using Reddy’s model also yields lower values for the natural frequencies. Table 3 lists the natural frequencies obtained using the full Bernoulli model and the second-order Reddy’s model ( - variant). Table 3. Natural frequencies of a heterogeneous beam ( = -1) obtained using the full Bernoulli model and the Reddy model Mode Bernoulli L/h = 5 Reddy L/h = 5 |w|max (Reddy L/h = 5 ) Bernoulli L/h = 10 Reddy L/h = 10 Bernoulli L/h = 20 Reddy L/h = 20 1 4.6143 4.4792 0.015 4.6484 4.6120 4.6570 4.6478 2 22.5405 19.7516 0.073 23.4913 22.5459 23.7304 23.4780 3 32.0069 31.8974 7.4 60.7498 55.7234 62.3378 60.7801 4 55.6523 44.2526 0.08 63.8790 63.7661 119.2573 114.0833 5 81.7920 72.3770 0.13 114.3008 99.3542 - - 6 98.5168 80.4830 4.9 - - - - S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. Figures 4 and 5 demonstrate the convergence of the developed finite element algorithm. Figure 4 shows how the obtained natural frequency values vary depending on the degree of discretization. Even with a fairly coarse finite element mesh, the accuracy of the solution is adequate. For example, dividing the beam into 8 nodal elements yields an accuracy of 0.03% for the first frequency and 0.4% for the third. The solutions obtained using the two variants of the algorithm ( - and -variants) differ to some extent, especially for the coarse mesh. However, in the -variant of the algorithm, a solution with the desired accuracy is achieved on a coarser mesh. Figure 5 shows that the results of the two algorithm variants converge as the number of nodal elements increases. The difference in the natural frequencies obtained using the - and -variants of the algorithm ( t and g ) is more significant for higher frequencies. The third mode is predominantly longitudinal (Table 3); the contribution of shear strain is negligible, and therefore the solutions obtained using the two variants of the algorithm are virtually identical. The results presented in Table 4 allow to assess the significance of accounting for the heterogeneity of beam properties when calculating natural frequencies. The table summarizes the solutions obtained for a beam that is homogeneous along its length ( = 0) with L/h 5 and for a beam whose properties vary in the longitudinal and transverse directions ( = -1). The solutions were obtained using the full Bernoulli model and the Reddy model. The table shows that accounting for the heterogeneous distribution of the physical and mechanical characteristics of a cantilever beam in the longitudinal direction leads to an increase in the natural frequencies and, consequently, to an increase in dynamic rigidity in the subresonance region. а b c Figure 4. Convergence of the numerical solution for (a) the first, (b) the third and (c) the fifth mode S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. Mode Bernoulli = -1 Reddy = -1 Bernoulli = 0 Reddy = 0 1 4.6143 4.4792 3.4330 3.3510 2 22.5405 19.7516 20.4895 17.9510 3 32.0069 31.8974 26.4676 26.3736 4 55.6523 44.2526 53.8033 42.7835 5 81.7920 72.3770 79.6095 71.1021 Table 4. Comparison of natural frequencies of a bidirectional ( = -1) and a longitudinally homogenous beam ( = 0) Figure 5. Relative frequency difference obtained using two variants of the algorithm. Digits indicate the mode number S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. S o u r c e: made by I.N. Vasserman, I.O. Glot, I.N. Shardakov. To test this finite-element implementation of the Reddy model, an alternative distribution of material properties presented in [28] was also considered, in which the Young’s modulus and mass density vary exponentially in both directions: E x, z Em exp px Lx 12 pz hz 12 ; x, z m exp px Lx 12 pz hz 12 . (46) The resulting dimensionless natural frequencies L2 m , (47) Table 5. Comparison of natural h Em frequencies with [28] Mode Karamanli [28] Present model 1 0.7055 0.7055 2 4.7302 4.7273 3 6.3414 6.3415 were compared with those calculated in [28] for the following parameter values: Em 210 GPa, m 7850 kg/m3, px 1, pz 1. The results are presented in Table 5. The deviation for the first three frequencies does not S o u r c e: made by I.N. Vasserman, exceed 0.1%. I.O. Glot, I.N. Shardakov. 4. Conclusion A finite-element algorithm is proposed for solving problems involving the static and dynamic deformation of bidirectional functionally graded beams and structures composed of them in terms of the Bernoulli and Reddy beam theories. As an example, the results of a numerical solution to the problem of natural modes and frequencies of a cantilever bidirectional functionally graded beam are presented. Based on these results, the following conclusions can be drawn: 1. The convergence of the obtained natural frequency values with an increase in the number of finite elements in the numerical solution has been established. In particular, for the simplified Bernoulli model, which does not account for the rotational inertia of the cross-sections or the cross-coupling of the longitudinal force and the bending moment, the obtained numerical values of the natural frequencies agree well with those of the analytical solution. 2. For short beams (L/h ≤ 5), accounting for rotational inertia leads to a significant reduction in natural frequencies, especially for the higher modes. 3. Numerical results that demonstrate the relationship between the natural frequencies and the asymmetry of the distribution of the physical properties of beam material in the cross-section relative to its centerline are presented. An interesting result is the obtained interdependence of the longitudinal and transverse vibrations of the beam when the stiffness properties are distributed asymmetrically relative to the centerline in the cross-sections. This interaction is unrelated to geometric or physical nonlinearity and should be distinguished from the effect of longitudinal loading on the natural frequency of transverse vibrations, which is closely associated with Euler instability. This result fundamentally distinguishes the vibrations of functionally graded beams from those of homogeneous beams. 4. Using Reddy’s higher-order model also results in a decrease in the calculated natural frequencies. For the lower frequencies of long beams, this difference is less noticeable. 5. It has been shown that an increase in the material porosity toward the free end of a cantilever beam leads to an increase in the natural frequencies and, consequently, to a decrease in dynamic displacements in the sub-resonance region. 6. Thus, the gradient spatial distribution of the beam material properties in two directions substantially broadens the range of control over the parameters of its dynamic and static deformation behavior.About the authors
Igor N. Vasserman
Institute of Continuous Media Mechanics, Ural Branch of the Russian Academy of Sciences (UB RAS)
Author for correspondence.
Email: igorw@icmm.ru
ORCID iD: 0000-0001-7799-7521
SPIN-code: 1961-6243
Candidate of Physical and Mathematical Sciences, Researcher
13A Lenin St, Perm, 614000, Russian FederationIrina O. Glot
Institute of Continuous Media Mechanics, Ural Branch of the Russian Academy of Sciences (UB RAS)
Email: glot@icmm.ru
ORCID iD: 0000-0002-2842-7511
SPIN-code: 3231-7996
Candidate of Physical and Mathematical Sciences, Senior Researcher
13A Lenin St, Perm, 614000, Russian FederationIgor N. Shardakov
Institute of Continuous Media Mechanics, Ural Branch of the Russian Academy of Sciences (UB RAS)
Email: shardakov@icmm.ru
ORCID iD: 0000-0001-8673-642X
SPIN-code: 4809-4776
Doctor of Physical and Mathematical Sciences, Chief Researcher
13A Lenin St, Perm, 614000, Russian FederationReferences
- Matveenko VP, Fedorov AYu, Shardakov IN. Analysis of stress singularities at singular points of elastic solids made of functionally graded materials. Doklady Physics. 2016;61(1):24–28. https://doi.org/10.1134/S1028335816010031 EDN: WQFCBN
- Fedorov AYu, Matveenko VP. Investigation of stress behavior in the vicinity of singular points of elastic bodies made of functionally graded materials. ASME Journal of Applied Mechanics. 2018;85(6):061008. https://doi.org/10.1115/1.4039619 EDN: OSQKPB
- Meguid SA, Wang XD, Jiang LY. On the dynamic propagation of a finite crack in functionally graded materials. Engineering Fracture Mechanics. 2002;69(14–16):1753–1768. https://doi.org/10.1016/S0013-7944(02)00046-2 EDN: AYEQFZ
- Bever MB, Duwez PE. Gradients in composite materials. Materials Science and Engineering. 1972;10:1–8. https://doi.org/10.1016/0025-5416(72)90059-6
- Miyamoto Y, Koizumi M, Yamada O. High-pressure self-combustion sintering for ceramics. Journal of the American Ceramic Society. 1984;67(11):224–225, https://doi.org/10.1111/j.1151-2916.1984.tb19488.x
- Niino M, Hirai T, Watanabe R. The functionally gradient materials, Journal of the Japan Society for Composite Materials. 1987;13(1):257–264. https://doi.org/10.6089/jscm.13.257
- BirmanV, Bird LW. Modelling and analysis of functionally graded materials and structures. ASME Applied Mechanics Reviews. 2007;60:195–216. https://doi.org/10.1115/1.2777164 EDN: SFXRPV
- Pompe W, Worch H, Epple M, Friess W, Gelinsky M, Greil P, Hempel U, Scharnweber D, Schulte K, Functionally graded materials for biomedical applications; Materials Science and Engineering: A; 2003;362(12):40-60. https://doi.org/
- 1016/S0921-5093(03)00580-X EDN: KLHYOD
- Oh SY, Librescu L, Song O. Vibration of turbomachinery rotating blades made-up of functionally graded materials and operating in a high temperature field. Acta Mechanica. 2003;166:69–87. https://doi.org/10.1007/s00707-003-0049-y EDN: EUNGWR
- Zhang N, Lu A, Charlie C, Li CC, Zhou J, Zhang X, Wang S, Xuguang Chen X. Support performance of functionally graded concrete lining. Construction and Building Materials. 2017;147:35-47. https://doi.org/10.1016/j.conbuildmat.2017.04.161
- Rubio WM, Vatanabe SL, Paulino GH, Silva ECN. Functionally graded piezoelectric material systems — A multi-physics perspective. In: Advanced Computational Materials Modeling: From Classical to Multi-Scale Techniques. 2010. p. 301–339. https://doi.org/10.1002/9783527632312.ch8
- Le CI, Dang ND, Bui TTH, Nguyen DK. Size-dependent large deformation analysis of microscale functionally graded beams and frames. Acta Mechanica. 2025;236:1991–2010. https://doi.org/10.1007/s00707-025-04265-1 EDN: CAGRMM
- Guo L, Zhao S, Guo Z, Yang J, Kitipornchai S. Functionally graded engineering structures: An overview. Engineering Structures. 2026;351:121995. https://doi.org/10.1016/j.engstruct.2025.121995 EDN: EJHCQB
- Torelli G, Fernández MG, Lees JM. Functionally graded concrete: Design objectives, production techniques and analysis methods for layered and continuously graded elements. Construction and Building Materials. 2020;242:118040. https://doi.org/10.1016/j.conbuildmat.2020.118040 EDN: HUGIZI
- Strieder E, Hilber R, Stierschneider E, Bergmeister K. FE-Study on the Effect of Gradient Concrete on Early Constraint and Crack Risk. Applied Sciences. 2018;8(2):246. https://doi.org/10.3390/app8020246
- Nian Y, Wan S, Avcar M, Yue R, Li M. 3D printing functionally graded metamaterial structure: Design, fabrication, reinforcement, optimization. International Journal of Mechanical Sciences. 2023;258:108580. https://doi.org/10.1016/
- j.ijmecsci.2023.108580 EDN: OXRZKO
- Ma Z, Liu W, Li W, Liu H, Song J. et al. Additive manufacturing of functional gradient materials: A review of research progress and challenges. Journal of Alloys and Compounds. 2024;971:172642. https://doi.org/10.1016/j.jallcom.
- 172642 EDN: GZIBML
- Lim S, Buswell RA, Le TT, Austin SA, Gibb AGF, Thorpe T. Developments in construction-scale additive manufacturing processes. Automation in Construction. 2012;21:262–268. https://doi.org/10.1016/j.autcon.2011.06.010
- Gosselin C, Duballet R, Roux P, Gaudillière N, Dirrenberger J, Morel P. Large-scale 3D printing of ultra-high performance concrete — a new processing route for architects and builders. Materials & Design. 2016;100:102–109. https://doi.org/10.1016/j.matdes.2016.03.097
- Tay YWD, Lim JH, Li M, Tan MJ. Creating functionally graded concrete materials with varying 3D printing para-meters. Virtual and Physical Prototyping. 2022;17(3):662–681. https://doi.org/10.1080/17452759.2022.2048521 EDN: GMGOJI
- Chen D, Gao K, Yang J, Zhang L. Functionally graded porous structures: Analyses, performances, and applications — A Review. Thin-Walled Structures. 2023;191:111046. https://doi.org/10.1016/j.tws.2023.111046 EDN: IMYZRY
- Shetty RA, Dushyanthkumar GL, Deepak SA. Classical and refined beam and plate theories: A brief technical review. Sahyadri International Journal of Research. 2018;4(2):27-32. Available from: https://sijr.sahyadri.edu.in/Articles/
- Volume%204/Issue%202/december_2018_6.pdf (accessed: 11.03.2026).
- Thai H-T, Vo TP. Bending and free vibration of functionally graded beams using various higher-order shear deformation beam theories. International Journal of Mechanical Sciences. 2012;62(1):57-66. https://doi.org/10.1016/j.ijmecsci.
- 05.014
- Şimşek M. Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories. Nuclear Engineering and Design. 2010;240(4):697-705. https://doi.org/10.1016/j.nucengdes.2009.12.013
- Emam SA. Analysis of shear-deformable composite beams in postbuckling. Composite Structures. 2011;94(1):
- -30. https://doi.org/10.1016/j.compstruct.2011.07.024
- Heyliger PR, Reddy JN. A Higher Order Beam Element for Bending and Vibration Problems. Journal of Sound and Vibration. 1988;126:309-326. https://doi.org/10.1016/0022-460X(88)90244-1
- Şimşek M. Bi-directional functionally graded materials (BDFGMS) for free and forced vibration of Timoshenko beams with various boundary conditions. Composite Structures. 2015;133:968–976. https://doi.org/10.1016/j.compstruct.2015.08.021
- Karamanli A. Free vibration analysis of two directional functionally graded beams using a third order shear deformation theory. Composite Structures. 2018;189:127–136. https://doi.org/10.1016/j.compstruct.2018.01.060
- Tang Y, Lv X, Yang T. Bi-directional functionally graded beams: asymmetric modes and nonlinear free vibration. Composites Part B: Engineering. 2019;156:319–331. https://doi.org/10.1016/j.compositesb.2018.08.140 EDN: YKHWNF
- Wang Z-H, Wang X-H, Xu G-D, Cheng S, Zeng T. Free vibration of two-directional functionally graded beams. Composite Structures. 2016;135:191-198. https://doi.org/10.1016/j.compstruct.2015.09.013
- Boggarapu V, Gujjala R, Ojha Sh, Acharya Sk, Venkateswara babu P, Chowdary S, kumar Gara D. State of the art in functionally graded materials. Composite Structures. 2021;262:113596. https://doi.org/10.1016/j.compstruct.2021.113596 EDN: YLQYWC
- Alshorbagy AE, Eltaher MA, Mahmoud FF. Free vibration characteristics of a functionally graded beam by finite element method. Applied Mathematical Modelling. 2011;35(1):412-425. https://doi.org/10.1016/j.apm.2010.07.006
- Pei YL, Geng PS, Li LX. A modified higher-order theory for FG beams. European Journal of Mechanics — A/Solids. 2018;72:186-197. https://doi.org/10.1016/j.euromechsol.2018.05.008
- Ferreira AJM. MATLAB codes for finite element analys. Springer; 2009. https://doi.org/10.1007/978-3-030-47952-7
Supplementary files










