Static and dynamic analyses of two-phase/multi-phase carbon nanotube-reinforced functionally graded composite beams via warping-included mixed finite element method

Abstract

This study enhances the application of cross-sectional warping considered mixed finite element (W-MFE) formulation to accurately determine natural vibration, static displacement response, and shear and normal stress evaluation with very close to the precision of solid finite elements (FEs) in two-phase/multi-phase functionally graded (FG) laminated composite beams strength using carbon nanotubes (CNTs). The principles of three dimensional (3D) elasticity theory are used to derive constitutive equations. The mixed finite element (MFE) method is improved by accounting for warping effects by displacement-based FEs within the cross-sectional domain. The MFE with two nodes has a total of 24 degrees of freedom. The two-phase material consists of a polymer matrix reinforced with aligned CNTs that are FG throughout the beam thickness. The multi-phase FG beam is modeled as a three-component composite material, consisting of CNTs, a polymer matrix, and fibers. The polymer matrix is reinforced by longitudinally aligned fibers and randomly dispersed CNT particles. The fiber volume fractions are considered to change gradually through the thickness of the beam following a power-law variation. The W-MFE achieves satisfactory results with fewer degrees of freedom than 3D solid FEs. Benchmark examples examine the effects of ply orientation, configuration, and fiber gradation on FG beam behavior.

Graphical abstract

Keywords

stress analysis / natural vibration analysis / CNT-reinforced composites / mixed finite element method / functionally graded material

1 Introduction

Polymer composites are extensively used in manufacturing and various industrial applications due to their low cost, ease of maintenance, and increased strength. However, these materials have significant drawbacks, such as delamination under high loads or extreme thermal conditions [1]. Functionally Graded Materials (FGMs) address these issues by gradually changing material properties in a specific direction. This gradation can be achieved by combining two or more parent materials using volume fractions or chemically treating a single material to alter its initial properties [2]. By gradually varying the composition of their constituents, FGMs can effectively mitigate stress concentrations, thermal stresses, and crack initiation at interfaces, improving the performance of traditional composites [3–5]. As a result, FGMs have found extensive applications in high-technology sectors, including defense, aerospace, biomedical, nuclear, energy conversion, electronics, and optics. Accurately predicting the behavior of structures made from advanced materials relies on mathematical modeling, with Partial Differential Equations (PDEs) serving as essential tools for understanding their response in natural and engineered systems. However, analytical solutions to PDEs are rarely feasible for complex, real-world problems, necessitating numerical methods. The finite element method (FEM) has long been a dominant approach in engineering applications, as it relies on discretizing the problem domain and approximating the solution with specifically chosen basis functions [6]. FEM remains a powerful and widely used method, particularly for solving complex engineering problems. Nevertheless, alternative approaches have been explored to address specific challenges, such as computational efficiency for high-dimensional problems or complex geometries. Among these, deep learning-based methods have demonstrated significant potential due to their ability to approximate complex functions and patterns flexibly and efficiently. In particular, Deep Neural Networks (DNNs) offer complementary capabilities to existing numerical methods by leveraging advancements in algorithmic efficiency and data-driven modeling. Samaniego et al. [7] emphasized the utility of the energetic formulation of PDEs in addressing mechanical problems, noting that the energy of a mechanical system serves as a natural loss function for machine learning methods to model such problems effectively. For instance, the deep collocation method proposed by Guo et al. [8] demonstrates the use of feedforward DNNs for approximating PDE solutions in thin plate bending problems. This method takes advantage of computational graphs and backpropagation algorithms inherent in deep learning. Similarly, Zhuang et al. [9] introduced a deep autoencoder-based energy method (DAEM) for analyzing the bending, vibration, and buckling of Kirchhoff plates. DAEM integrates the higher-order continuity of autoencoders with the minimum total potential principle, offering an unsupervised learning framework that provides advantages over conventional methods.

Carbon nanotubes (CNTs) [10,11] or fibers [12,13] are commonly employed as reinforcement phases in polymer composites. CNTs are renowned for their excellent thermo-electro-mechanical properties, making them ideal reinforcements for composite materials [14]. Surface treatments or functionalization of CNTs enhance interfacial bonding with the matrix, thereby improving load transfer efficiency [15–17]. This improved load transfer significantly enhances the overall mechanical properties of the composites. Similarly, fibers significantly improve the tensile strength, compressive strength, and stiffness of composite materials. The orientation, type, and volume fraction of fibers are crucial in defining these strength characteristics, and their gradation through the thickness is important for optimizing performance [18–21]. To further understand the factors influencing the mechanical properties of composites, it is important to consider the distribution and alignment of reinforcing materials. Mechanical performance in nanocomposites depends on these factors, as imperfections like agglomeration or waviness can reduce stiffness and load-bearing capacity, thereby diminishing overall mechanical properties. Rafiee and collaborators comprehensively investigated the effects of CNT agglomeration and waviness on the material properties of carbon nanotube-reinforced composites (CNTRCs) [22–25] along with extensive research, including multiscale modeling and stochastic approaches [26–30]. For nanoclay/polymer composites [31] and graphene/polymer composites [32], similar challenges, such as particle dispersion and alignment, are likely to influence their material properties, highlighting the broader implications of these factors on the mechanical properties of various nanocomposite systems.

Given the increasing reliance on FGMs, it is essential to have a thorough understanding and analyzing the mechanical behavior and performance of functionally graded (FG) beams. Accurate prediction of the behavior of FG beams is crucial to prevent structural failure, as they are employed in structures subjected to diverse loading conditions, including transverse, vibration, and in-plane loads. The orientation of the CNTs within the polymer matrix significantly influences the mechanical characteristics of CNTRCs. The literature identifies two primary states of CNT orientation in two-phase composites: aligned CNTs [33–35] and randomly oriented CNTs [36,37]. It is evident that most research on FG beams reinforced by CNTs predominantly focuses on aligned or randomly oriented CNT configurations. Several shear deformation theories have been presented in the literature to perform static or dynamic analyses of beam-type structures. Among the theories are those such as Classical Beam Theory (CBT), First-Order Shear Deformation Theory (FSDT), Higher-Order Shear Deformation Theory (HSDT) [38–40], and Zigzag Theories (ZZT) [41–43]. Some of the studies related to two-phase composite beams consisting of aligned or randomly oriented CNTs are as follows: CBT, also referred to as Euler-Bernoulli Beam Theory, is based on the assumption that plane sections remain planar and orthogonal to the mid-surface during deformation, which makes it highly suitable for the analysis of slender beams. Heshmati and Yas [44] performed the free vibration characteristics of FG-randomly oriented CNTRC beams by FEM, based on the CBT. The FSDT, widely known as Timoshenko Beam Theory, was introduced to address the inadequacy of CBT. FSDT enhances CBT by accounting for rotational inertia and first-order shear effects. Yas and Heshmati [45] investigated the vibration characteristics of FG-randomly oriented CNTRC beams under the action of a moving load employing FEM. FSDT and CBT were used to examine the dynamic behavior of the beam. Ansari et al. [46] analyzed the forced vibration behavior of nanocomposite beams reinforced with single-walled aligned CNTs, employing the FSDT and accounting for nonlinear geometric effects based on von Kármán theory. Jam and Kiani [47] investigated the low-velocity impact behavior of FG-aligned CNTRC beams based on the FSDT in a thermal environment. Using the generalized differential quadrature technique, Kamarian et al. [48] studied the free vibration analysis of FG-randomly oriented CNTRC beams based on the FSDT, resting on a Pasternak foundation. Ebrahimi and Farazamandnia [49] used the differential transform method to examine the free vibration characteristics of FG-aligned CNTRC sandwich beams based on FSDT in a thermal environment. Borjalilou et al. [50] studied the bending, buckling, and free vibration behaviors of FG-aligned CNTRC nanobeams based on the FSDT and non-local elasticity theory, considering small-scale effects. Talebi et al. [51] used the Chebyshev-Ritz method to examine the thermal free vibration behavior of sandwich piezoelectric agglomerated beams based on FSDT strengthened by randomly oriented CNTs, considering the pyroelectric effect. HSDT were introduced as an advancement over FSDT to improve modeling precision and more accurately represent shear deformations in beams and plates, eliminating the requirement for a shear correction coefficient. HSDT can capture transverse deformations through thickness, accounting for cross-sectional warping by incorporating higher-order shear terms expressed using polynomial, trigonometric, exponential functions, etc. Lin and Xiang [52] analyzed the free vibration behavior of FG-aligned CNTRC beams based on first- and third-order beam theories using the p-Ritz method. Jedari Salami [53] presented an extended high-order sandwich panel theory to examine the bending behavior of FG-aligned CNTRC sandwich beams with a soft core. Biswas and Datta [54] investigated the free vibration analysis of FG-aligned CNTRC beams using FEM based on refined shear deformation theories. Madenci [55] developed a mixed FEM to investigate the free vibration characteristics of FG-aligned CNTRC beams based on trigonometric shear deformation theory (SDT). Belarbi et al. [56] analyzed the bending and buckling behavior of FG-aligned CNTRC beams with the FEM, employing a hyperbolic SDT. ZZT use piecewise continuous functions to refine the in-plane displacement field and accurately represent cross-sectional deformation caused by nonuniform transverse shear stress in multi-layered composites. They ensure the continuity of the shear stress field across layers while keeping the number of kinematic variables constant. Chalak et al. [57] studied the free vibration and modal stress analysis of FG-aligned CNTRC beam by using FEM based on higher order ZZT. A survey of the existing literature shows that hybrid composite materials are widely used to enhance composite structures by incorporating multi-phase materials in distinct forms, as highlighted in several studies [58–67]. The following hybrid/multiscale composite beam studies are noted in the recent literature. He et al. [68] proposed an analytical formulation for the large-amplitude free and forced vibration response of CNT/fiber/polymer laminated multiscale composite beams based on CBT and von Kármán geometric nonlinearity. Afshin and Yas [69] investigated the dynamic and buckling analysis of a polymer hybrid composite beam with variable thickness based on FSDT using the differential quadrature method. Alambeigi et al. [70] conducted a comprehensive study on the free and forced vibration analysis of a sandwich beam with a porous core and shape memory alloy hybrid composite face layers on Vlasov’s foundation. The governing equations of motion were derived using Hamilton’s principle and FSDT. An analytical approach was employed to solve these equations, utilizing Navier’s method. Using a refined beam theory, Dabbagh et al. [71] investigated the thermal buckling analysis of agglomerated multiscale hybrid nanocomposites. Daikh et al. [72] analyzed the static bending response of laminated randomly oriented FG-CNT/fiber-reinforced composite beams resting on an elastic foundation. Li et al. [73] investigated the nonlinear vibrational characteristics of multiscale composite beams resting on a foundation with nonlinear softening behavior using Reddy’s third-order SDT, which considers nonlinear geometric effects based on von Kármán theory.

The existing literature on composite beams reinforced with CNTs primarily concentrates on two-phase composite materials, with limited investigation of multi-phase composite materials. This highlights the need for further research into the static and dynamic behaviors of CNT-reinforced FG layered beams, particularly within a unified framework that incorporates shear-locking-free and warping-considered mixed finite element (FE) method (W-MFEM). The W-MFEM, which exhibits convergence behavior comparable to 3D FEs, utilizes a two-node beam element with 24 degrees of freedom (DOF). This formulation is derived from Timoshenko beam theory, enhanced by a condensed three dimensional (3D) elasticity framework. The method accounts for cross-sectional warping through a displacement-based FE formulation integrated within the MFE approach. Normal stresses are evaluated using constitutive relations incorporating warping effects, while shear stresses are derived from the warping-integrated deformation functions and axial rotations. Recent studies have validated the accuracy and computational efficiency of W-MFEM: Aribas et al. [74] analyzed the behavior of static loads and normal/shear stress distributions of curved beams with super-elliptical geometry having axially FGM using W-MFEM, highlighting its accuracy in capturing complex stress distributions. Aribas et al. [75] further emphasized the significance of cross-sectional warping and demonstrated the enhanced accuracy of W-MFEM compared to FSDT in analyzing the dynamic behavior of planar curved beams reinforced with transversely FGM. Ermis [76] also studied the dynamic response of spatially curved beams reinforced with FG graphene platelets subjected to oscillatory excitation forces. This study extends the application of W-MFEM to two-phase and multi-phase CNT-reinforced FG layered beams, addressing the aforementioned research gaps. The following sections detail the methodology and results of this investigation. Section 2 presents the micromechanical modeling of two-phase and three-phase CNT-reinforced composite beams, including the homogenization methods that predict their effective material properties. Section 3 outlines the mixed FE formulation and its application to static and dynamic analyses of laminated FG composite beams. In Section 4, numerical examples validate the proposed method through comparisons with existing literatures [77–81] and 3D FE analyses. Additionally, parametric studies investigate the effects of material gradation, beam configurations, and ply orientations on the free vibration, static, and stress behaviors of the CNT/polymer/fiber FG laminated composite beams.

2 The micromechanical model for the mechanical properties of the two-phase/multi-phase composite beam

2.1 Two-phase composite beam

The two-phase FG-CNTRC beams, extensively investigated in the literature by various researchers, including Refs. [77–81], are considered. This two-phase material is characterized by a polymer matrix enriched solely with aligned CNTs ([33,34], Fig.1). Classical linear-graded patterns of aligned CNT, including UD, X, O, and V, are used to design these beams (Fig.1). In these configurations, the aligned CNT reinforcing phase exhibits a linear variation throughout the thickness of the beam, representing the graded structures described in previous studies as follows:

(1)\{\begin{array}{rl}& \mathrm{U}\mathrm{D}:V_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(\mathit{\text{z}}\right)={\tilde{V}}_{\mathrm{C}\mathrm{N}\mathrm{T}},\\ & \mathrm{F}\mathrm{G}-\mathrm{X}:V_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(\mathit{\text{z}}\right)={\tilde{V}}_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(\frac{4\left|\mathit{\text{z}}\right|}{h}\right),\\ & \mathrm{F}\mathrm{G}-\mathrm{O}:V_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(\mathit{\text{z}}\right)={\tilde{V}}_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(2-\frac{4\left|\mathit{\text{z}}\right|}{h}\right),\\ & \mathrm{F}\mathrm{G}-\mathrm{V}:V_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(\mathit{\text{z}}\right)={\tilde{V}}_{\mathrm{C}\mathrm{N}\mathrm{T}}\left(1+\frac{2\mathit{\text{z}}}{h}\right),\end{array}

where {\tilde{V}}_{\mathrm{C}\mathrm{N}\mathrm{T}} is the volume fraction of aligned CNT. The effective material properties (E_{11}, E_{22}= E_{33}, G_{12}= G_{13}, and {\upsilon }_{12}={\upsilon }_{13}) of the CNTRC beams are calculated by the rule of mixture in the reference papers [77,79,80] as follows,

(2)E_{11}= {\eta }_1{V}_{\mathrm{C}\mathrm{N}\mathrm{T}}{E}_{11}^{\mathrm{C}\mathrm{N}\mathrm{T}}+V_{\mathrm{M}}{E}_{\mathrm{M}},

(3)\{\begin{array}{rl}& \frac{{\eta }_2}{E_{22}}=\frac{V_{\mathrm{C}\mathrm{N}\mathrm{T}}}{E_{22}^{\mathrm{C}\mathrm{N}\mathrm{T}}}+\frac{V_{\mathrm{M}}}{E_{\mathrm{M}}},\\ & \frac{{\eta }_3}{G_{12}}=\frac{V_{\mathrm{C}\mathrm{N}\mathrm{T}}}{G_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}}+\frac{V_{\mathrm{M}}}{G_{\mathrm{M}}},\end{array}

(4){\upsilon }_{12}=V_{\mathrm{C}\mathrm{N}\mathrm{T}}{\upsilon }_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}+V_{\mathrm{M}}{E}_{\mathrm{M}},

(5)\rho =V_{\mathrm{C}\mathrm{N}\mathrm{T}}{\rho }_{\mathrm{C}\mathrm{N}\mathrm{T}}+V_{\mathrm{M}}{\rho }_{\mathrm{M}},

where {\eta }_i(i =1,2,3) indicates the aligned CNT efficiency parameter.

In this study, the effective shear modulus is considered G_{12}=G_{23} as stated in Karamanli and Vo [79]. The effective Poisson’s ratio is considered as {\upsilon }_{23}= \left(E_{22}/E_{11}\right){\upsilon }_{12} in Ref. [15]. The material properties of the matrix and CNTs used in this study for the two-phase material are as follows [79]: the material properties of aligned CNT is E_{11}^{\mathrm{C}\mathrm{N}\mathrm{T}}=600\mathrm{G}\mathrm{P}\mathrm{a}, E_{22}^{\mathrm{C}\mathrm{N}\mathrm{T}}=10\mathrm{G}\mathrm{P}\mathrm{a}, G_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}=17.2\mathrm{G}\mathrm{P}\mathrm{a}, {\upsilon }_{12}^{\mathrm{C}\mathrm{N}\mathrm{T}}=0.19, {\rho }_{\mathrm{C}\mathrm{N}\mathrm{T}}=1400 \text{kg/}{\text{m}}^3. The following material properties characterize the matrix material: E_{\mathrm{M}}=2.5\mathrm{G}\mathrm{P}\mathrm{a}, {\upsilon }_{\mathrm{M}}= 0.3, {\rho }_{\mathrm{M}}=1190\text{kg/}{\text{m}}^3.

2.2 Multi-phase composite beam

Advanced multiscale composite materials are extensively used in numerous industrial applications, such as turbine blades, cutting tools, and aircraft engines. The beams under consideration are made of three-phase composite layers, consisting of FG-CNT/polymer/fiber composite materials. These materials comprise a polymer matrix strengthened with straight-aligned glass fibers and randomly distributed CNTs (Fig.2). A multiscale method is required to evaluate the engineering constants that govern the effective material properties of the multiphase composite. The evaluation is conducted in a two-step process presented in Ref. [64]. Initially, Eshelby-Mori-Tanaka (EMT) approach [82] is utilized at the nanoscale to determine the characteristics of the polymer matrix enhanced with only CNTs. These CNTs are described as particles with transverse isotropy that are uniformly distributed within the matrix and oriented randomly [83–85]. Then, the Hahn method is employed to integrate these properties and the characteristics of the reinforcing fibers [86].

Following the EMT approach, the evaluation of the mechanical properties of the CNT-strengthened polymer matrix are evaluated through bulk K_{\mathrm{M}}^{\ast } and shear moduli G_{\mathrm{M}}^{\ast }, as detailed below

(6)\{\begin{array}{rl}& K_{\mathrm{M}}^{\ast }=K_{\mathrm{M}}+\frac{V_{\mathrm{C}}\left({\delta }_{\mathrm{C}}-3K_{\mathrm{M}}{\alpha }_{\mathrm{C}}\right)}{3\left(V_{\mathrm{M}}+V_{\mathrm{C}}{\alpha }_{\mathrm{C}}\right)},\\ & G_{\mathrm{M}}^{\ast }=G_{\mathrm{M}}+\frac{V_{\mathrm{C}}\left({\eta }_{\mathrm{C}}-2G_{\mathrm{M}}{\beta }_{\mathrm{C}}\right)}{2\left(V_{\mathrm{M}}+V_{\mathrm{C}}{\beta }_{\mathrm{C}}\right)},\end{array}

where K_{\mathrm{M}}= E_{\mathrm{M}}/3\left(1-2{\upsilon }_{\mathrm{M}}\right) and G_{\mathrm{M}}= E_{\mathrm{M}}/2\left(1+{\upsilon }_{\mathrm{M}}\right), being E_{\mathrm{M}}, {\upsilon }_{\mathrm{M}} the Young’s modulus and the Poisson’s ratio of the pure polymer matrix. The volume fraction of the CNTs V_{\mathrm{C}} is defined by the following expression, which is a function of the densities of the two constituents {\rho }_{\mathrm{C}}, {\rho }_{\mathrm{M}} and the mass fraction of the reinforcing phase w_{\mathrm{C}}:

(7)V_{\mathrm{C}}= {\left(\frac{{\rho }_{\mathrm{C}}}{w_{\mathrm{C}}{\rho }_M}-\frac{{\rho }_{\mathrm{C}}}{{\rho }_{\mathrm{M}}}+1\right)}^{-1},

being V_{\mathrm{M}}= 1-V_{\mathrm{C}} the volume fraction of the matrix. To address the necessary definitions Eq. (6), The subsequent parameters are necessary:

(8){\alpha }_{\mathrm{C}}= \frac{3\left(K_{\mathrm{M}}+G_{\mathrm{M}}\right)+k_{\mathrm{C}}-l_{\mathrm{C}}}{3\left(G_{\mathrm{M}}+k_{\mathrm{C}}\right)},

(9)\begin{array}{rl}{\beta }_{\mathrm{C}}=& \frac{1}{5}(\frac{4G_{\mathrm{M}}+2k_{\mathrm{C}}+l_{\mathrm{C}}}{3\left(G_{\mathrm{M}}+k_{\mathrm{C}}\right)}+\frac{4G_{\mathrm{M}}}{G_{\mathrm{M}}+p_{\mathrm{C}}}\\ & +\frac{2\left(G_{\mathrm{M}}\left(6K_{\mathrm{M}}+8G_{\mathrm{M}}\right)\right)}{G_{\mathrm{M}}\left(3K_{\mathrm{M}}+G_{\mathrm{M}}\right)+m_{\mathrm{C}}\left(3K_{\mathrm{M}}+7G_{\mathrm{M}}\right)}),\end{array}

(10){\delta }_{\mathrm{C}}= \frac{1}{3}\left(n_{\mathrm{C}}+2l_{\mathrm{C}}+\frac{\left(2k_{\mathrm{C}}+l_{\mathrm{C}}\right)\left(3K_{\mathrm{M}}+2G_{\mathrm{M}}-l_{\mathrm{C}}\right)}{G_{\mathrm{M}}+k_{\mathrm{C}}}\right),

(11)\begin{array}{rl}{\eta }_{\mathrm{C}}=& \frac{1}{5}(\frac{2}{3}\left(n_{\mathrm{C}}-l_{\mathrm{C}}\right)+\frac{8G_{\mathrm{M}}{p}_{\mathrm{C}}}{G_{\mathrm{M}}+p_{\mathrm{C}}}+\frac{2\left(k_{\mathrm{C}}-l_{\mathrm{C}}\right)\left(2G_{\mathrm{M}}+l_{\mathrm{C}}\right)}{3\left(G_{\mathrm{M}}+k_{\mathrm{C}}\right)}\\ & +\frac{8m_{\mathrm{C}}{G}_{\mathrm{M}}\left(3K_{\mathrm{M}}+4G_{\mathrm{M}}\right)}{3K_{\mathrm{M}}\left(m_{\mathrm{C}}+G_{\mathrm{M}}\right)+G_{\mathrm{M}}\left(7m_{\mathrm{C}}+G_{\mathrm{M}}\right)}),\end{array}

where k_{\mathrm{C}}, l_{\mathrm{C}}, m_{\mathrm{C}}, n_{\mathrm{C}}, p_{\mathrm{C}} define the Hill’s elastic moduli of CNT. Subsequently, the Youngs’s modulus E_{\mathrm{M}}^{\ast }, the density {\rho }_{\mathrm{M}}^{\ast }, and the Poisson’s ratio {\upsilon }_{\mathrm{M}}^{\ast } of the enriched matrix are defined as

(12)\{\begin{array}{rl}& E_{\mathrm{M}}^{\ast }=\frac{9K_{\mathrm{M}}^{\ast }{G}_{\mathrm{M}}^{\ast }}{3K_{\mathrm{M}}^{\ast }+G_{\mathrm{M}}^{\ast }},\\ & {\upsilon }_{\mathrm{M}}^{\ast }=\frac{3K_{\mathrm{M}}^{\ast }-2G_{\mathrm{M}}^{\ast }}{6K_{\mathrm{M}}^{\ast }+2G_{\mathrm{M}}^{\ast }},\\ & {\rho }_{\mathrm{M}}^{\ast }=\left({\rho }_{\mathrm{C}}-{\rho }_{\mathrm{M}}\right)V_{\mathrm{C}}+{\rho }_{\mathrm{M}},\end{array}

where V_{\mathrm{F}}= 1-V_{\mathrm{M}}^{\ast }. In this context, the composite maintains its isotropic properties, a result of the random distribution of the nanoparticles, as demonstrated in the study by Shi et al. [87]. At this point, the incorporation of straight reinforcing fibers into the composite can be achieved using Hahn’s homogenization procedure [86]. Initially, it is necessary to compute the fiber volume fraction {\tilde{V}}_{\mathrm{F}} in relation to the corresponding density {\rho }_{\mathrm{F}} and mass fraction w_{\mathrm{F}} as described by the following equation:

(13){\tilde{V}}_{\mathrm{F}}={\left(\frac{{\rho }_{\mathrm{F}}}{w_{\mathrm{F}}{\rho }_M^{\ast }}-\frac{{\rho }_{\mathrm{F}}}{{\rho }_{\mathrm{M}}^{\ast }}+1\right)}^{-1}.

In this study, the longitudinally aligned straight reinforcing fibers, having a non-uniform distribution, are considered to be located through the beam’s thickness coordinate z using the relation,

(14)V_{\mathrm{F}}= {\tilde{V}}_{\mathrm{F}}{f}_i^{\left(k\right)}\left(\mathit{\text{z}}\right),

(15)f_1^{\left(k\right)}\left(\mathit{\text{z}}\right)={\left(\frac{\mathit{\text{z}}-{\mathit{\text{z}}}_k}{{\mathit{\text{z}}}_{k+1}-{\mathit{\text{z}}}_k}\right)}^{\alpha }\mathrm{o}\mathrm{r}{f}_2^{\left(k\right)}\left(\mathit{\text{z}}\right)={\left(\frac{{\mathit{\text{z}}}_{k+1}-\mathit{\text{z}}}{{\mathit{\text{z}}}_{k+1}-{\mathit{\text{z}}}_k}\right)}^{\alpha },

for \alpha \in {\mathbb{R}}^{+}. According to the proposed power-law functions in Eq. (15), a variety of complex patterns can be achieved by changing the exponent \alpha of the power-law functions f_i^{\left(k\right)}\left(\mathit{\text{z}}\right), for i=1 and 2 which are the generic layers of composite beams. The power-law functions are shown in Fig.3 for the different exponent values of \alpha =0.1,0.3,0.5,1, 2,5, and 20. The blue arrows in Fig.3 indicate the direction of the increase in exponent \alpha.

For the first power law function f_1^{\left(k\right)}\left(\mathit{\text{z}}\right), at the bottom of the generic layer (when \mathit{\text{z}}={\mathit{\text{z}}}_k) there are no fibers and at the top of the generic layer (when \mathit{\text{z}}={\mathit{\text{z}}}_{k+1}) the bottom of the generic layer has the fiber volume fraction is equal {\tilde{V}}_F (Fig.3(a), Eq. (14)). For the second power law function f_2^{\left(k\right)}\left(\mathit{\text{z}}\right), at the bottom of the generic layer (when \mathit{\text{z}}={\mathit{\text{z}}}_k) has the same fiber volume fraction is equal {\tilde{V}}_{\mathrm{F}} and at the top of the generic layer (when \mathit{\text{z}}={\mathit{\text{z}}}_{k+1}) there are no fibers (Fig.3(b)). As the value of the exponent \alphaincreases, for f_1^{\left(k\right)}\left(\mathit{\text{z}}\right) the top of the generic region becomes fiber-rich and for f_2^{\left(k\right)}\left(\mathit{\text{z}}\right) the bottom of the generic region becomes fiber-rich.

The kth three-phase layer has engineering constants expressed as

(16)E_{11}^{\left(k\right)}=E_{11}^{\mathrm{F}}{V}_{\mathrm{F}}+ E_{\mathrm{M}}^{\ast }{V}_{\mathrm{M}}^{\ast },

(17){\upsilon }_{12}^{\left(k\right)}={\upsilon }_{12}^{\mathrm{F}}{V}_{\mathrm{F}}+ {\upsilon }_{\mathrm{M}}^{\ast }{V}_{\mathrm{M}}^{\ast },

(18)G_{12}^{\left(k\right)}=G_{13}^{\left(k\right)}=\frac{V_{\mathrm{F}}+{\mathrm{\Delta }}_1{V}_{\mathrm{M}}^{\ast }}{V_{\mathrm{F}}/G_{12}^F+{\mathrm{\Delta }}_1{V}_{\mathrm{M}}^{\ast }/G_{\mathrm{M}}^{\ast }} ,

(19)G_{23}^{\left(k\right)}=\frac{V_{\mathrm{F}}+{\mathrm{\Delta }}_2{V}_{\mathrm{M}}^{\ast }}{V_{\mathrm{F}}/G_{23}^{\mathrm{F}}+{\mathrm{\Delta }}_2{V}_{\mathrm{M}}^{\ast }/G_{\mathrm{M}}^{\ast }} ,

(20)E_{22}^{\left(k\right)}=\frac{4E_{11}^{\left(k\right)}{K}_{\mathrm{T}}^{\left(k\right)}{G}_{23}^{\left(k\right)}}{E_{11}^{\left(k\right)}{K}_{\mathrm{T}}^{\left(k\right)}+G_{23}^{\left(k\right)}\left(E_{11}^{\left(k\right)}+4K_{\mathrm{T}}^{\left(k\right)}{\left({\upsilon }_{12}^{\left(k\right)}\right)}^2\right)},

where the following parameters are needed:

(21)\{\begin{array}{rl}& K_{\mathrm{T}}={\displaystyle \frac{V_{\mathrm{F}}+{\mathrm{\Delta }}_{\mathrm{K}}{V}_{\mathrm{M}}^{\ast }}{{\displaystyle \frac{V_{\mathrm{F}}}{K_{\mathrm{F}}}}+{\mathrm{\Delta }}_{\mathrm{K}}{\displaystyle \frac{V_{\mathrm{M}}^{\ast }}{K_{\mathrm{M}}^{\ast }}}}},\\ & {\mathrm{\Delta }}_{\mathrm{K}}={\displaystyle \frac{1+{\displaystyle \frac{G_{\mathrm{M}}^{\ast }}{K_{\mathrm{F}}}}}{2\left(1-{\upsilon }_{\mathrm{M}}^{\ast }\right)}},\\ & {\mathrm{\Delta }}_1={\displaystyle \frac{1+{\displaystyle \frac{G_{\mathrm{M}}^{\ast }}{G_{12}^{\mathrm{F}}}}}{2}},\\ & {\mathrm{\Delta }}_2={\displaystyle \frac{3-4{\upsilon }_M^{\ast }+{\displaystyle \frac{G_{\mathrm{M}}^{\ast }}{G_{23}^{\mathrm{F}}}}}{4\left(1-{\upsilon }_{\mathrm{M}}^{\ast }\right)}},\end{array}

(22)E_{33}^{\left(k\right)}=E_{22}^{\left(k\right)},{\upsilon }_{13}^{\left(k\right)}={\upsilon }_{12}^{\left(k\right)},

(23){\upsilon }_{23}^{\left(k\right)}={\upsilon }_{12}^{\mathrm{F}}{V}_{\mathrm{F}}+ {\upsilon }_{\mathrm{M}}^{\ast }{V}_{\mathrm{M}}^{\ast }\left(\frac{1+{\upsilon }_{\mathrm{M}}^{\ast }-{\upsilon }_{12}^{\left(k\right)}{E}_{\mathrm{M}}^{\ast }/E_{11}^{\left(k\right)}}{1+{\left({\upsilon }_{\mathrm{M}}^{\ast }\right)}^2+{\upsilon }_{\mathrm{M}}^{\ast }{\upsilon }_{12}^{\left(k\right)}{E}_{\mathrm{M}}^{\ast }/E_{11}^{\left(k\right)}}\right).

In the proposed approach, the reinforcing fibers’ properties are described concerning of Young’s moduli E_{11}^{\mathrm{F}}, E_{22}^{\mathrm{F}}, shear moduli G_{12}^{\mathrm{F}}, G_{23}^{\mathrm{F}}, Poisson’s ratio {\upsilon }_{12}^{\mathrm{F}}, {\upsilon }_{23}^F, and bulk modulus K_{\mathrm{F}}= {\displaystyle \frac{E_{22}^{\mathrm{F}}}{3\left(1-2{\upsilon }_{23}^{\mathrm{F}}\right)}}. If the reinforcing fibers are characterized by isotropic features, these properties can be simplified as follows:

(24)\{\begin{array}{rl}& E_{\mathrm{F}}=E_{11}^{\mathrm{F}}=E_{22}^{\mathrm{F}},\\ & G_{\mathrm{F}}=G_{12}^{\mathrm{F}}=G_{23}^{\mathrm{F}},\\ & {\upsilon }_{\mathrm{F}}={\upsilon }_{12}^{\mathrm{F}}={\upsilon }_{23}^{\mathrm{F}}.\end{array}

It should be kept in mind that the parameters mentioned in Eqs. (16)–(23) vary with the thickness due to specified fiber volume fraction. The material properties of the straight isotropic fibers and the polymer matrix are given in Tab.1, while the characteristics of single-walled CNT with 20 as chiral index are provided in Tab.2 based on of Hill’s elastic moduli. The overall material density of the composite layer is determined using the rule of mixtures. For the subsequent applications, the required mass fractions w_{\mathrm{F}}=0.85 and w_{\mathrm{C}}=0.05 specified in Eqs. (7)–(13) are used (unless stated otherwise).

3 Warping considered mixed finite element formulation for multi-layer functionally graded composite beam

3.1 The constitutive equations

According to the generalized Hooke’s law, the relationship between the stresses \mathit{\sigma } and strains \mathit{\epsilon } in an elastic continuum is defined by the three-dimensional stress-strain components as \mathit{\sigma }=\mathit{E} \mathit{\epsilon } where \mathit{E} is the elasticity matrix. An axis transformation on the elasticity matrix can be obtained by a matrix operation \overline{\mathit{E}}={\mathit{T}}^{-1}\mathit{E}{\mathit{T}}^{-\mathrm{T}} where \mathit{T} is the axis transformation matrix, and the superscript ‘T’ denotes the matrix transpose operator [88]. Finally, letting x, y, and z be the Cartesian Coordinates and the beam axis be x, the reduction of the three-dimensional elasticity theory to the beam theory is achieved by letting {\sigma }_y={\sigma }_{\mathit{\text{z}}}={\tau }_{y\mathit{\text{z}}}=0 [89], the constitutive equation in matrix form yields

(25){\left\{\begin{array}{l}{\sigma }_x\\ {\tau }_{\mathit{\text{z}}x}\\ {\tau }_{xy}\end{array}\right\}}_k={\left[\begin{array}{ccc}{\overline{\beta }}_{11}& 0& {\overline{\beta }}_{13}\\ 0& {\overline{\beta }}_{22}& 0\\ {\overline{\beta }}_{13}& 0& {\overline{\beta }}_{33}\end{array}\right]}_k{\left\{\begin{array}{l}{\epsilon }_x\\ {\gamma }_{\mathit{\text{z}}x}\\ {\gamma }_{xy}\end{array}\right\}}_k= {\overline{\mathit{\beta }}}_k{\left\{\begin{array}{c}{\epsilon }_x\\ {\gamma }_{\mathit{\text{z}}x}\\ {\gamma }_{xy}\end{array}\right\}}_k,

where {\overline{\mathit{\beta }}}_k is the transformed reduced elasticity matrix [90] (Appendix A in Supplementary materials), and k is the number of each laminae. The strain components at any point on the beam cross-section can be expressed in terms of the displacements \mathit{u}\left(u_x,u_y,u_{\mathit{\text{z}}}\right) and corresponding cross-sectional rotations \mathbf{\Omega }\left({\Omega }_x,{\Omega }_y,{\Omega }_{\mathit{\text{z}}}\right) along the beam axis, based on the Timoshenko beam theory, as follows [91]:

(26)\{\begin{array}{l}{\epsilon }_x=u_{x,x}^{\ast },\\ {\gamma }_{xy}=u_{x,y}^{\ast }+u_{y,x}^{\ast }=u_{y,x}-\mathit{\text{z}}{\Omega }_{x,x},\\ {\gamma }_{\mathit{\text{z}}x}=u_{x,\mathit{\text{z}}}^{\ast }+u_{\mathit{\text{z}},x}^{\ast }=u_{\mathit{\text{z}},x}+y{\Omega }_{x,x}.\end{array}

Single-layer constitutive equations using kinematic relations yield [92],

(27){\left\{\begin{array}{l}{\sigma }_x\\ {\tau }_{\mathit{\text{z}}x}\\ {\tau }_{xy}\end{array}\right\}}_k={\overline{\mathit{\beta }}}_k{\left\{\left\{\begin{array}{l}{u}_{x,x}\\ u_{\mathit{\text{z}},x}\\ u_{y,x}\end{array}\right\}+\mathit{\text{z}}\left\{\begin{array}{c}{\Omega }_{y,x}\\ 0\\ -{\Omega }_{x,x}\end{array}\right\}+y\left\{\begin{array}{c}-{\Omega }_{\mathit{\text{z}},x}\\ {\Omega }_{x,x}\\ 0\end{array}\right\}\right\}}_k,k=1,\dots ,N_{\mathrm{L}}.

Commas in subscripts denote partial derivatives, and the total number of laminae is represented by N_{\mathrm{L}}. The nodal forces and moments are derived by an analytical integration of stresses for individual layer across the cross-sectional thickness, as detailed in Refs. [90,93].

(28)\{\begin{array}{rl}& F_x=\sum _{k=1}^{N_{\mathrm{L}}}\left[{\int }_{A_k}{\left({\sigma }_x\right)}_k\mathrm{d}{A}_k\right],\\ & F_y=\sum _{k=1}^{N_{\mathrm{L}}}\left[{\int }_{A_k}{\left({\tau }_{xy}\right)}_k\mathrm{d}{A}_k\right],\\ & F_{\mathit{\text{z}}}=\sum _{k=1}^{N_{\mathrm{L}}}\left[{\int }_{A_k}{\left({\tau }_{\mathit{\text{z}}x}\right)}_k\mathrm{d}{A}_k\right],\\ & M_x=\sum _{k=1}^{N_{\mathrm{L}}}\left[-{\int }_{A_k}\mathit{\text{z}}{\left({\tau }_{xy}\right)}_k\mathrm{d}{A}_k+{\int }_{A_k}y{\left({\tau }_{x\mathit{\text{z}}}\right)}_k\mathrm{d}{A}_k\right],\\ & M_y=\sum _{k=1}^{N_{\mathrm{L}}}\left[{\int }_{A_k}\mathit{\text{z}}{\left({\sigma }_x\right)}_k\mathrm{d}{A}_k\right],\\ & M_{\mathit{\text{z}}}=\sum _{k=1}^{N_{\mathrm{L}}}\left[-{\int }_{A_k}y{\left({\sigma }_x\right)}_k\mathrm{d}{A}_k\right],\end{array}

where F_x is axial force, F_y, and F_{\mathit{\text{z}}} are the shear forces, M_x represents the torque, M_y and M_{\mathit{\text{z}}} denote bending moments, respectively. A_k denotes the area of the cross-section for individual lamina. Thus, the constitutive equations in a matrix form yield,

(29)\begin{array}{rl}\left\{\begin{array}{l}{F}_x\\ F_y\\ F_{\mathit{\text{z}}}\\ M_x\\ M_y\\ M_{\mathit{\text{z}}}\end{array}\right\}=& \sum _{k=1}^{N_{\mathrm{L}}}{\left[\begin{array}{ll}{\mathit{E}}^{\mathrm{m}}& {\mathit{E}}^{\mathrm{m}\mathrm{f}}\\ {\mathit{E}}^{\mathrm{f}\mathrm{m}}& {\mathit{E}}^{\mathrm{f}}\end{array}\right]}_k\left[\begin{array}{l}{u}_{x,x}\\ u_{y,x}\\ u_{\mathit{\text{z}},x}\\ {\Omega }_{x,x}\\ {\Omega }_{y,x}\\ {\Omega }_{\mathit{\text{z}},x}\end{array}\right\}\\ =& \sum _{k=1}^{N_{\mathrm{L}}}{\left[\begin{array}{cccccc}{E}_{11}^{\mathrm{m}}& E_{12}^{\mathrm{m}}& 0& E_{11}^{\mathrm{m}\mathrm{f}}& E_{12}^{\mathrm{m}\mathrm{f}}& 0\\ E_{21}^{\mathrm{m}}& E_{22}^{\mathrm{m}}& 0& E_{21}^{\mathrm{m}\mathrm{f}}& E_{22}^{\mathrm{m}\mathrm{f}}& 0\\ 0& 0& E_{33}^{\mathrm{m}}& 0& 0& 0\\ E_{11}^{\mathrm{f}\mathrm{m}}& E_{12}^{\mathrm{f}\mathrm{m}}& 0& E_{11}^{\mathrm{f}}& E_{12}^{\mathrm{f}}& 0\\ E_{21}^{\mathrm{f}\mathrm{m}}& E_{22}^{\mathrm{f}\mathrm{m}}& 0& E_{21}^{\mathrm{f}}& E_{22}^{\mathrm{f}}& 0\\ 0& 0& 0& 0& 0& E_{33}^{\mathrm{f}}\end{array}\right]}_k\left\{\begin{array}{c}{u}_{x,x}\\ u_{y,x}\\ u_{\mathit{\text{z}},x}\\ {\Omega }_{x,x}\\ {\Omega }_{y,x}\\ {\Omega }_{\mathit{\text{z}},x}\end{array}\right\},\end{array}

where {\mathit{E}}^{\mathrm{m}}, {\mathit{E}}^{\mathrm{m}\mathrm{f}}, {\mathit{E}}^{\mathrm{f}\mathrm{m}}, and {\mathit{E}}^{\mathrm{f}} are obtained from Eq. (28). {\mathit{E}}^{\mathrm{m}} and {\mathit{E}}^{\mathrm{f}} are elastic stiffness matrices related to force and moments, respectively. {\mathit{E}}^{\mathrm{m}\mathrm{f}} and {\mathit{E}}^{\mathrm{f}\mathrm{m}} (transpose of {\mathit{E}}^{\mathrm{m}\mathrm{f}}) represent the coupling elastic stiffness matrices. The torsional stiffness of composite cross-sections with non-circular geometry, accounting for warping effects, is provided in Ref. [94],

(30){\stackrel{︵}{GI}}_{\mathrm{t}}={\int }_A\mathit{d}\mathit{G}(\mathit{d}-\mathit{B} \widehat{\text{ψ}})\mathrm{d}A,

in this framework, the coordinate vector for the mesh on the cross-section is referred to as \mathit{d}={\{\begin{array}{ll}\mathit{\text{z}}& -y\end{array}\}}^{\mathrm{T}} [93,95]. I_{\mathrm{t}} is the torque moment of inertia, the matrix \mathit{G} represents the shear modulus of the mesh element. \widehat{\text{ψ}} corresponds to the warping magnitude vector corresponding to the nodes, and the partial derivatives vector of the shape functions are denoted as \mathit{B}. The constitutive equation for a beam yields,

(31)\left\{\begin{array}{l}{u}_{x,x}\\ u_{y,x}\\ u_{\mathit{\text{z}},x}\\ {\mathrm{\Omega }}_{x,x}\\ {\mathrm{\Omega }}_{y,x}\\ {\mathrm{\Omega }}_{\mathit{\text{z}},x}\end{array}\right\}=\left[\begin{array}{ll}{\mathit{C}}^{\mathrm{m}}& {\mathit{C}}^{\mathrm{m}\mathrm{f}}\\ {\mathit{C}}^{\mathrm{f}\mathrm{m}}& {\mathit{C}}^{\mathrm{f}}\end{array}\right]\left\{\begin{array}{c}{F}_x\\ F_y\\ F_{\mathit{\text{z}}}\\ M_x\\ M_y\\ M_{\mathit{\text{z}}}\end{array}\right\} ,

where the matrix components of the compliance matrix \mathit{C}={\mathit{E}}^{-1} are {\mathit{C}}^{\mathrm{m}}, {\mathit{C}}^{\mathrm{f}}, {\mathit{C}}^{\mathrm{m}\mathrm{f}}, and {\mathit{C}}^{\mathrm{f}\mathrm{m}} (transpose of {\mathit{C}}^{\mathrm{m}\mathrm{f}}), respectively.

3.2 Field equations and functionals

The field equations for spatially curved isotropic Timoshenko beams [96] have been adapted to develop the governing equations for laminated composite straight beams,

(32)\begin{array}{rl}& \begin{array}{r}-{\mathit{F}}_{,x}-\mathit{q}+\left({\displaystyle \sum _{k=1}^{N_{\mathrm{L}}}{\rho }_k{A}_k}\right)\ddot{\mathit{u}}=\mathbf{0}\\ -{\mathit{M}}_{,x}-\mathbf{i}\times \mathit{F}-\mathit{m}+\left({\displaystyle \sum _{k=1}^{N_{\mathrm{k}}}{\rho }_k{\mathit{I}}_k}\right)\ddot{\mathbf{\Omega }}=\mathbf{0}\end{array}\}\\ & \text{equations of motion},\end{array}

(33)\begin{array}{r}{\mathit{u}}_{,x}+\mathbf{i}\times \mathbf{\Omega }-{\mathit{C}}^{\mathrm{m}}\mathit{F}-{\mathit{C}}^{\mathrm{m}\mathrm{f}}\mathit{M}=\mathbf{0}\\ {\mathbf{\Omega }}_{,x}-{\mathit{C}}^{\mathrm{f}\mathrm{m}}\mathit{F}-{\mathit{C}}^{\mathrm{f}}\mathit{M}=\mathbf{0}\end{array}\}\text{constitutive equations},

(34)\begin{array}{rr}-\mathit{u}+\hat{\mathbf{u}}=\mathbf{0};& -\mathbf{\Omega }+\hat{\mathbf{\Omega }}=\mathbf{0}\\ \mathit{F}-\hat{\mathit{F}}=\mathbf{0};& \mathit{M}-\hat{\mathit{M}}=\mathbf{0}\end{array}\}\text{boundary conditions}.

In these equations, \mathit{u}\left(u_x,u_y,u_{\mathit{\text{z}}}\right) denotes the displacement vector, \mathbf{\Omega }\left({\Omega }_x,{\Omega }_y,{\Omega }_{\mathit{\text{z}}}\right) represents the rotation vector, \mathit{F}\left(F_x,F_y,F_{\mathit{\text{z}}}\right) shows the force vector, \mathit{M}\left(M_x,M_y,M_{\mathit{\text{z}}}\right) denotes the moment vector. The total number of laminae is N_{\mathrm{L}}, and k identifies the number of each lamina. The overall material density of composite kth laminae is represented by {\rho }_k, with the cross sectional area A_k, and the moment of inertia {\mathit{I}}_k. The vectors \mathit{q} and \mathit{m} indicate the distributed external loads in the form of forces and moments, respectively. The acceleration terms are as follows: \ddot{\mathit{u}} and \ddot{\mathbf{\Omega }}. The given values at the boundary are denoted with hats. Using the Gateaux differential and potential operator [97,98], the necessary form of the functional for the static and free vibration analysis is derived based on Eqs. (32)–(34) as follows.

Considering the static analysis, it is clear that \ddot{\mathit{u}}=\ddot{\mathbf{\Omega }}=0, the functional for static analysis is obtained as follows:

(35)\begin{array}{rl}{\mathit{I}}_{\mathrm{s}}\left(\mathit{y}\right)=& -\left[\mathit{u},{\mathit{F}}_{,x}\right]+[\mathbf{i}\times \mathbf{\Omega },\mathit{F}]-\left[{\mathit{M}}_{,x},\mathbf{\Omega }\right]-\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{m}}\right)\mathit{F},\mathit{F}\right]\\ & -\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{m}\mathrm{f}}\right)\mathit{M},\mathit{F}\right]-\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{f}\mathrm{m}}\right)\mathit{M},\mathit{F}\right]\\ & -\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{f}}\right)\mathit{M},\mathit{M}\right]-[\mathit{q},\mathit{u}]-[\mathit{m},\mathbf{\Omega }]+{⟨(\mathit{F}-\hat{\mathit{F}}),\mathit{u}⟩}_{\sigma }\\ & +{⟨(\mathit{M}-\hat{\mathit{M}}),\mathbf{\Omega }⟩}_{\sigma }+{⟨\mathit{F},\hat{\mathit{u}}⟩}_{\epsilon }+{⟨\mathit{M},\hat{\mathbf{\Omega }}⟩}_{\epsilon }.\end{array}

Considering free vibration analysis under harmonic motion, it follows that, the external force/moment vectors are zero. The acceleration components can be expressed as [\ddot{\mathit{u}},\ddot{\mathit{u}} ]=- {\omega }^2[\mathit{u},\mathit{u}] and [\ddot{\mathbf{\Omega }},\ddot{\mathbf{\Omega }} ]=- {\omega }^2[\mathbf{\Omega },\mathbf{\Omega }], where \omega represents the natural circular frequency. The functional for free vibration analysis is obtained as follows:

\begin{array}{rl}{\mathit{I}}_{\omega }\left(\mathit{y}\right)=& -\left[\mathit{u},{\mathit{F}}_{,x}\right]+[\mathbf{i}\times \mathbf{\Omega },\mathit{F}]-\left[{\mathit{M}}_{,x},\mathbf{\Omega }\right]-\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{m}}\right)\mathit{F},\mathit{F}\right]\\ & -\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{m}\mathrm{f}}\right)\mathit{M},\mathit{F}\right]-\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{f}\mathrm{m}}\right)\mathit{M},\mathit{F}\right]\\ & -\frac{1}{2}\left[\left({\mathit{C}}^{\mathrm{f}}\right)\mathit{M},\mathit{M}\right]-\frac{1}{2}\left(\sum _{k=1}^{N_{\mathrm{L}}}{\rho }_k{A}_k\right){\omega }^2[\mathit{u},\mathit{u}]\end{array}

(36)\begin{array}{rl}& -\frac{1}{2}{\omega }^2\left[\left(\sum _{k=1}^{N_{\mathrm{L}}}{\rho }_k{\mathit{I}}_k\right)\mathbf{\Omega },\mathbf{\Omega }\right]+{⟨(\mathit{F}-\hat{\mathit{F}}),\mathit{u}⟩}_{\sigma }\\ & +{⟨(\mathit{M}-\hat{\mathit{M}}),\mathbf{\Omega }⟩}_{\sigma }+{⟨\mathit{F},\hat{\mathit{u}}⟩}_{\epsilon }+{⟨\mathit{M},\hat{\mathbf{\Omega }}⟩}_{\epsilon }.\end{array}

The inner product is defined through square bracket notation. The subscript ‘\epsilon’ denote geometric boundary condition, while the subscript \sigma corresponds to dynamic boundary conditions, respectively.

3.3 Mixed finite element formulation

The two-noded straight beam mixed FE element is constructed by two linear shape functions {\phi }_i\left(\xi \right), {\phi }_j\left(\xi \right). where \xi \in \left[0,1\right] is the local coordinate, and the subscripts i and j are the node numbers of the element. The nodal unknow vectors are \mathit{u}\left(u_x,u_y,u_{\mathit{\text{z}}}\right), \mathit{\Omega }\left({\Omega }_x,{\Omega }_y,{\Omega }_{\mathit{\text{z}}}\right), \mathit{F}\left(F_x,F_y,F_{\mathit{\text{z}}}\right), \mathit{M}\left(M_x,M_y,M_{\mathit{\text{z}}}\right). Finally, the mixed FE matrix {\mathit{k}}^{\mathrm{e}} and mass matrix {\mathit{m}}^{\mathrm{e}} are obtained as follows:

(37){\mathit{k}}^{\mathrm{e}}=\left[\begin{array}{cccc}\mathbf{0}& \mathbf{0}& {\mathit{k}}_{\mathrm{u}\mathrm{F}}& \mathbf{0}\\ & \mathbf{0}& {\mathit{k}}_{\mathrm{\Omega }\mathrm{F}}& {\mathit{k}}_{\mathrm{\Omega }\mathrm{M}}\\ & & {\mathit{k}}_{\mathrm{F}\mathrm{F}}& {\mathit{k}}_{\mathrm{F}\mathrm{M}}\\ & \text{sym.}& & {\mathit{k}}_{\mathrm{M}\mathrm{M}}\end{array}\right],{\mathit{m}}^{\mathrm{e}}= \left[\begin{array}{cccc}{\mathit{m}}_{\mathrm{u}\mathrm{u}}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ & {\mathit{m}}_{\mathrm{\Omega }\mathrm{\Omega }}& \mathbf{0}& \mathbf{0}\\ & & \mathbf{0}& \mathbf{0}\\ & \text{sym.}& & \mathbf{0}\end{array}\right].

Submatrices of Eq. (37) exists in Appendix B in Supplementary materials. Detailed mixed FE formulations for static and free vibration analyses of beams, plates, and shells exists in Refs. [98–101].

3.4 Calculation of natural frequencies

Determining the natural frequencies of vibration in a structural system is formulated as a standard eigenvalue problem, expressed as \left|\left(\mathit{K}-{\omega }^2\mathit{M}\right)\right|=0, where and \omega is the natural angular frequency of the system. By using the conventional FE assembly technique, including the necessary boundary conditions, the system matrix K and the system mass matrix \mathit{M} can be constructed. Consequently, the problem is reduced to an eigenvalue problem:

(38)\left(\left[\begin{array}{ll}{\mathit{K}}_{11}& {\mathit{K}}_{12}\\ {\mathit{K}}_{12}^{\mathrm{T}}& {\mathit{K}}_{22}\end{array}\right]-{\omega }^2\left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathit{M}\end{array}\right]\right)\left\{\begin{array}{l}\mathit{R}\\ \mathit{D}\end{array}\right\}=\left\{\begin{array}{l}\mathbf{0}\\ \mathbf{0}\end{array}\right\},

where \mathit{D} is the displacements plus rotations column vector, and \mathit{R} is the stress resultants (forces plus moments) column vector. For the free vibration analysis by mixed FE formulation, elimination of \mathit{R} in Eq. (38) reduces to a standard eigenvalue problem [101].

(39)\left({\mathit{K}}^{\ast }-{\omega }^2\mathit{M}\right) \left\{\mathit{D}\right\}=\left\{\mathbf{0}\right\} ,

where {\mathit{K}}^{\ast }={\mathit{K}}_{22}-{\mathit{K}}_{12}^{\mathrm{T}}{\mathit{K}}_{11}^{-1}{\mathit{K}}_{12} is the condensed system matrix. All numerical computations were carried out using a custom-built computer program based on the mixed FE algorithm described above, implemented in the FORTRAN programming language.

3.5 Stress analysis

The calculation of normal stress distribution across the cross-section is derived from the constitutive equations in conjunction with the relevant curvatures in Eq. (27), while shear stress distribution is determined through the both axial rotations and warping function of the MFEs framework, as outlined in Ref. [94],

(40)\mathit{\tau }={\Omega }_{x,x}\mathit{G}(\mathrm{\nabla }\psi -\mathit{d}),

where the stress vector is represented as \mathit{\tau }={\left\{\begin{array}{ll}{\tau }_{xy}& {\tau }_{xz}\end{array}\right\}}^{\mathrm{T}}, the shear modulus of mesh elements is denoted as \mathit{G}, the vector of nodal coordinates in the mesh is shown as \mathit{d}. The axial rotation corresponds to {\Omega }_{x,x}. \mathrm{\nabla } denotes the in-plane differential operator, and \psi represents the function of warping [93]. The flowchart of the MFE algorithm is given in Fig.4.

4 Numerical examples

4.1 Convergence and comparison

4.1.1 Functionally graded reinforcing phase (two-phase material)

The purpose of this example is to determine the natural frequencies and static characteristics of FG configurations with a satisfactory level of precision, showcasing the capability of the W-MFEs. To achieve this, the organization of this example is as follows. First, a FE mesh and a number of layer convergence analyses are performed to assess the convergence of the W-MFE formulation for free vibration/static behaviors of FG-CNTRC beams with various distribution patterns, including UD, X, O, and V. Secondly, the W-MFEM convergence performances and SOLID186 elements (3D quadratic solid FEs, ANSYS software [102]) are investigated. Lastly, 3D quadratic solid FEs for free vibration/static analyses of the CNTRC beam, based on various CNT distribution patterns (namely, UD, O, X, and V), are compared with the results of W-MFE formulation and existing literature. The slenderness ratio of the beam L/h=15. The CNT volume fraction is {\tilde{V}}_{\mathrm{C}\mathrm{N}\mathrm{T}}=0.12 where the efficiency parameters are as follows: {\eta }_1= 1.2833, {\eta }_2={\eta }_3=1.0556. The boundary conditions are defined as clamped-clamped (CC) and clamped-free (CF). Detailed discussions are provided for free vibration and static analysis in the following subsections as follows:

1) Free vibration analysis

A FE mesh and the number of layer convergence analyses are performed to reveal the performance of the W-MFE formulation in reflecting the dynamic characteristics of FG CNTRC beam. The reference results of UD, X, O and V CNTRC beam are taken from Yas and Samadi [77], Vo-Duy et al. [78], Karamanli and Vo [79], and Garg et al. [80]. Yas and Samadi [77] and Vo-Duy et al. [78] applied the FSDT. Karamanli and Vo [79] adopted the shear and normal deformation theory (SNDT), and Garg et al. [80] used the higher-order ZZT. The non-dimensional out-of-plane circular natural frequency parameter is obtained by using the formula \tilde{\omega }=\omega L\sqrt{{\rho }_{\mathrm{M}}\left(1-{\upsilon }_{\mathrm{M}}^2\right)/E_{\mathrm{M}}}.

During the verification examples, a FE mesh convergence study is first conducted, along with an analysis of the number of elements, for the first three non-dimensional out-of-plane natural frequencies of the W-MFEs. For this purpose, the FG-V distribution is chosen. The analysis is carried out using 132 to 972 DOFs, corresponding to 10 to 80 elements, with the number of layers fixed at N_{\mathrm{L}}= 50 for CC and CF boundary conditions. Convergence results for the FG-V distribution frequencies indicate that 972 DOFs and the chosen number of layers are sufficient for the W-MFEs. Subsequently, a layer convergence analysis is performed for W-MFEs using 10 to 50 layers, with the necessary DOFs set to 972. The results of this analysis are tabulated in Tab.3. The percent differences are calculated based on the non-dimensional natural frequencies of the FG-V CNTRC beam for the W-MFEs (N_{\mathrm{L}}= 50, DOFs: 972) in comparison to those obtained using SOLID186 (N_{\mathrm{L}}= 50, DOFs: 2436000), as shown in Tab.3. The fundamental natural frequencies obtained from the W-MFEs are higher than those from the SOLID186 elements. As shown in Tab.3, the results for both the W-MFEs and the 3D quadratic solid FEs are in satisfactory consistency.

Secondly, the results of the present formulation (W-MFE) and 3D quadratic solid FEs for the first three natural frequencies of the CNTRC beam, based on different CNT distribution patterns (namely, UD, O, X, and V), are compared with the findings presented by Yas and Samadi [77], Vo-Duy et al. [78], Karamanli and Vo [79], and Garg et al. [80] in Tab.4–Tab.7, respectively. The results of WFEM in Tab.4–Tab.7 are represented N_L = 50 and DOFs: 972. When compared to the literature results for the clamped free (CF) boundary condition, the fundamental natural frequency results of the W-MFEs (Tab.4–Tab.7) are very closely agreed with those of the 3D elements (SOLID186) for all distribution patterns considered (UD, O, X and V). It can be seen that for the first three natural frequencies of CC boundary condition, the W-MFEM results converge to the SOLID186 results in terms of engineering precision. Specifically, the literature results for the first natural frequency show convergence from above compared to the SOLID186 results, while the W-MFEM results show convergence from below compared to the SOLID186 results for the UD and X-type distribution patterns (Tab.4 and Tab.5). For the O and V distribution patterns, both the MFEM and literature results exhibit convergence from above compared to the SOLID186 results (Tab.6 and Tab.7). As the dominant frequency is approached, the convergence of the W-MFEM results with SOLID186 becomes more pronounced for all-considered distribution patterns.

2) Static analysis

The FG-CNTRC beams are subjected to a uniform distributed load. The central deflection (x={\displaystyle \frac{L}{2}}) is non-dimensional as {\tilde{u}}_{\mathit{\text{z}}}={\displaystyle \frac{100E_{\mathrm{M}}{h}^{3}w}{q{L}^4}}{u}_{\mathit{\text{z}}}. The reference results are taken from the studies as follows: Kumar and Srinivas [81], and Karamanli and Vo [79] in Tab.8. Kumar and Srinivas applied FSDT. Karamanli and Vo adopted SNDT. The results of the non-central deflection of the CNTRC beam in Tab.8 belong to the warping-included mixed FEs (N_{\mathrm{L}}= 50, DOFs: 972) and SOLID186 elements (N_{\mathrm{L}}= 50, DOFs: 2436000). The analysis results from the W-MFEs and 3D solid FEs for central deflection of the CNTRC beam, based on various CNT distribution patterns (namely, UD, O, X, and V), are compared with the findings presented by Kumar and Srinivas [81] and Karamanli and Vo [79] in Tab.8. Compared with the literature results, the W-MFEs closely match those of the 3D solid FEs (SOLID186).

4.1.2 Functionally graded reinforcing phase (three-phase material)

This example evaluates the precision of the W-MFE model developed herein for analyzing the static and free vibration behavior of three-phase CNT/polymer/fiber FG composite beams. The linear distribution of straight reinforcing fibers (Subsection 2.2) is applied through the thickness of a three-phase composite FG beam with a two-layer configuration, where each layer has the same thickness. The first layer is composed of generic layer 1, while the second layer is composed of generic layer 2 (Fig.3). To achieve this, the natural frequency and static response results determined via W-MFEs are evaluated against those computed through 3D solid FEs (SOLID186) in ANSYS software. The geometrical properties of the straight beam are {\lambda }_{\mathrm{s}}= L/ h=20 and w/h=1 where w=0.04\mathrm{m}, with the beam clamped at both ends.

1) Free vibration analysis

A mesh convergence study is conducted for W-MFEs on the first three natural frequencies of a three-phase FG composite beam with a two-layer configuration, featuring either symmetric (0°/0°) or antisymmetric (0°/90°) layups, and a total of N_{\mathrm{L}}= 40. The analysis is performed using 252 to 972 DOFs, corresponding to 20 to 80 elements. Convergence results for the frequencies of both symmetric and antisymmetric layups indicate that 972 DOFs are sufficient for the W-MFEs. Subsequently, a layer convergence analysis is performed for W-MFEs using 20 to 80 layers, with the necessary DOFs set to 972. The results are tabulated in Tab.9, demonstrating that sufficient precision is achieved for the natural frequencies with N_{\mathrm{L}}= 80.

A mesh refinement analysis of the quadratic 3D solid FEs is conducted to examine the first three natural frequencies of the three-phase composite FG beam with respect to the total layer N_{\mathrm{L}}= 40. The DOFs are 390240, 1493280, and 3703692, respectively. The results of 3D solid FEs are tabulated in Tab.10. The percent differences (Tab.10) are calculated by normalizing the results obtained at 1493280 DOFs to those obtained at 3703692 DOFs. 0.03% is obtained as the absolute maximum percent difference. Tab.10 shows that the number of DOFs required for convergence is 3703692 for SOLID186. Next, a comparison analysis shows the convergence performance of the W-MFEs and 3D solid FEs with 40 layers in Tab.10. The percent differences are calculated by normalizing the frequencies obtained from W-MFEs (a total of 972 DOFs) to those obtained from SOLID 186 elements (DOFs: 3703692). 0.97% is obtained as the absolute maximum percent difference. The W-MFE results show a high level of consistency with the solutions from 3D solid FEs. In addition, the convergence of the W-MFE element is obtained with fewer DOFs than that of SOLID186.

2) Static analysis

The three-phase composite CNT/polymer/fiber FG beam with a two-layer configuration is under a uniformly distributed load having a load magnitude of q_{\mathit{\text{z}}} =10\mathrm{N} /\mathrm{m}. First, an FE mesh and a number of layer convergence analyses are performed for W-MFEs to obtain sufficiently convergent results of the displacement, support reactions, and stress components. The displacement is obtained at the midpoint on the axis of the beam. The stress components are calculated through the thickness of the beam at the axis x=0.1L. The FE convergence analysis is performed for DOFs 252, 492, 732, and 972 with a 160-layer number. The percent differences are calculated by normalizing the results obtained at 732 DOFs to those obtained at 972 DOFs. Calculated at 0.11% is the absolute maximum percent variation. It can be concluded that the number of DOFs required for convergence is 972 for W-MFEs. Next, the number of layers convergence analyses is performed for 40, 80, and 160 layers with 972 DOFs. The percent differences (Tab.11) are calculated by normalizing the results obtained at N_{\mathrm{L}}= 80 to the results obtained at N_{\mathrm{L}}= 160. The maximum absolute percent difference is 0.03% when considering the maximum normal/shear stress values. From Tab.11, it can be concluded that the number of N_{\mathrm{L}} required for convergence is 160 for W-MFEs.

Then, an FE mesh convergence analysis of SOLID186 is performed to obtain sufficiently convergent results of the displacement, support reactions, and stress components, as shown in Tab.12. The FE convergence analysis is conducted with DOFs 390240, 468720, 3518412 and 4286892 with a 40-layer number. The percent differences (Tab.12) are calculated by normalizing the results obtained at 3518412 DOFs with a 40-layer number to those obtained at 4286892 DOFs with the same number of layers. The maximum absolute percent difference is 0.004%. Tab.12 shows that the number of DOFs required for convergence is 4286892 DOFs with a 40-layer number for SOLID186. The maximum number of layers chosen for SOLID186 is 40 due to computer capacity (Intel(R) Core^TM i7-10750H CPU, 2.60GHz, 64GB RAM). Therefore, a comparison analysis is performed to show the convergence performance of the W-MFEs and quadratic 3D solid FEs with 40 layers. The percent differences are calculated by normalizing the results obtained from W-MFES (DOFs:972, Tab.11) to those obtained from SOLID 186 elements (DOFs: 4286892, Tab.12). The percent differences of F_{\mathit{\text{z}}}^{\mathrm{A}},M_x^{\mathrm{A}}, u_{\mathit{\text{z}}}^{\mathrm{M}}, {\sigma }_x(0.1,- 0.3),{\tau }_{x\mathit{\text{z}}} \left(0.1,0\right) and {\sigma }_x(0.5,- 0.3) are 0.00%, 0.19%, –1.11%, 0.36%, –5.81%, –0.56%, respectively. A detailed comparative analysis is performed by examining the absolute maximum shear stress values at the different cross-sections along the beam axis. For this purpose, the normalized length coordinates are chosen between the range 0.1\le \overline{x}\le 0.2. Fig.5 shows the absolute maximum shear stress through the thickness (\overline{\mathit{\text{z}}}). Normalization of WMFE results is performed with respect to those obtained from quadratic 3D solid FEs. The maximum percent difference is −5.81% shown in Fig.5. The ratio between the absolute maximum normal stress and shear stress of W-MFE (N_{\mathrm{L}}= 40) in Tab.11 is nearly 5.2. Additionally, it is noted that the shear stress values are lower than the normal stress values along the beam axis. These results show good consistency between the W-MFE and quadratic 3D solid FEs results. Based on the outcomes of the convergence analyses for both static and free vibration analyses, 80 warping-included MFEs are used to discretize the beams along the length.

4.2 Parametric solutions-free vibration

This example investigates how different lamination schemes and variations in fiber distribution (compositional gradient exponent) affect the natural frequencies of a three-phase CNT/polymer/fiber FG composite beam. In parametric analysis, two through-thickness configurations, Configuration 1 and Configuration 2 (Fig.6), are used, with each featuring two layers of equal thickness in either a symmetric (0°/0°) or antisymmetric (0°/90°) lamination scheme. The analysis also examines different span-to-thickness ratios ({\lambda }_{\mathrm{s}}= L/ h). The span ratio is set to {\lambda }_{\mathrm{s}}= 10, 15, and 20, successively where L is the length and h is the thickness. The beam length varies based on the h values (h=0.04\mathrm{m}), while the width is constant at w=0.04\mathrm{m}. The compositional gradient exponent \alpha is successively set to 0.1, 0.5, 1, 2, and 3 for each case. In Configuration 1, the mid-section of the cross-section has the highest glass fiber volume fraction, while the top and bottom sections contain no glass fiber (Fig.6(a)). In contrast, in Configuration 2, the top and bottom sections have the highest glass fiber volume fraction, with no glass fiber present in the mid-section (Fig.6(b)). Also, decreasing the compositional gradient exponent leads to an increased glass fiber volume fraction throughout the thickness of the beam for both configurations, resulting in stiffer FG material constituents through the thickness.

The first three non-dimensional transverse natural frequencies of a three-phase composite beam with a two-layer Configuration 1 and 2, featuring 0°/0° and 0°/90° layups for different span-to-thickness ({\lambda }_{\mathrm{s}}) are tabulated in Tab.13–Tab.15, respectively. The non-dimensional natural frequency is calculated by {\tilde{f}}_i=f_i{L}^2\sqrt{{\rho }_{\mathrm{M}}A/E_{\mathrm{M}}{I}_y}, where i=1,2, and 3.

The effect of the compositional gradient index (\alpha): Reduction of the composition gradient index increases the stiffness of the FG material through thickness (Fig.6) and results in higher natural frequencies in Configurations 1 and 2, both for the 0°/0° and 0°/90° lamination schemes. The corresponding numerical results are shown in Tab.13–Tab.15. The percentage difference in the natural frequencies is calculated by normalizing the results for α = 0.5, 1, 2, and 3 with respect to the results for α = 0.1 by using \mathrm{d}\mathrm{i}\mathrm{f}\mathrm{f}.\mathrm{\%} =\left(1-{\tilde{f}}_i^{\alpha =j}/{\tilde{f}}_i^{\alpha =0.1}\right)\times 100, where i=1,2, and 3, j=0.5,1,2, and 3. Fig.7(a) and Fig.7(b) show the percentage differences plotted against the compositional gradient index for Configuration 1, considering both 0°/0° and 0°/90° lamination schemes. Similarly, Fig.7(c) and Fig.7(d) show the same comparison for Configuration 2 using the corresponding lamination schemes. Also, the fundamental natural frequencies of Configuration 1 are more significantly influenced by variations in the compositional gradient index than Configuration 2. Specifically, Configuration 1, with the fiber-rich region centrally located within the cross-section, exhibits a more pronounced sensitivity to compositional gradient changes than Configuration 2, where the fiber-rich regions are situated at the upper and lower parts of the cross-section. When the higher stiffness regions (fiber-rich parts) are positioned at the upper and lower sections, as in Configuration 2, the structure becomes less sensitive to variations in the compositional gradient index than Configuration 1.

Additionally, for Configuration 1, the maximum reduction in natural frequencies occurs at the highest span-to-thickness ratio of 20, while the minimum reduction is observed at the lowest span-to-thickness ratio of 10, for both 0°/0° and 0°/90° lamination schemes across all compositional gradient indices. In contrast, for Configuration 2, the maximum increase in natural frequencies is observed at the lowest span-to-thickness ratio of 10, while the minimum increase occurs at the highest span-to-thickness ratio of 20.

The effect of the configurations: Configuration 2 consistently exhibits the highest values for the first three transverse natural frequencies (Tab.13–Tab.15). The percentage differences are plotted against the compositional gradient index in Fig.8(a) for the 0°/0° lamination scheme and in Fig.8(b) for the 0°/90° lamination scheme, respectively. The maximum reduction in natural frequencies is seen in the first mode, while the minimum reduction occurs in the third mode, with all other parameters held constant. Additionally, the reduction in natural frequencies becomes more substantial as the compositional gradient index increases in the 0°/0° lamination scheme compared to the 0°/90° lamination scheme.

4.3 Parametric solutions-static analysis

The influences of configuration and compositional gradient exponent on the vertical displacements and normal/shear stresses of three-phase CNT/polymer/fiber FG beams subjected to a uniformly distributed load q_{\mathit{\text{z}}} =10\mathrm{N} /\mathrm{m} are investigated. The beam is fixed at both ends. Two configurations are considered: Configuration 1 and Configuration 2, corresponding to 0°/0° and 0°/90° layups, respectively. The compositional gradient exponent of fiberglass is analyzed for values of 0.5, 1, and 2. Detailed analyses of the vertical displacement, normal stress, and shear stress are presented in Subsubsections 4.3.1, 4.3.2, and 4.3.3, respectively.

4.3.1 The vertical displacement

The vertical displacements u_{\mathit{\text{z}}} at the midspan of a three-phase composite beam having a two-layer Configuration 1 and 2 with 0°/0° and 0°/90° layups for different span-to-thickness ratio {\lambda }_{\mathrm{s}}= 10, 15, and 20 are tabulated in Tab.16. The non-dimensional maximum vertical displacement of the beam is obtained by {\tilde{u}}_{\mathit{\text{z}}}=\left(100E_{\mathrm{M}}{h}^{3}w/q{L}^4\right)u_{\mathit{\text{z}}}.

The effect of the compositional gradient index (α): a decrease in the compositional gradient index increased the stiffness of the FG beam and reduced the vertical displacements at the midspan of the beams (Tab.16). The percent differences are calculated by comparing the results for α = 1 and α = 0.5 with respect to α = 2 for each configuration and lamination scheme. Decreasing the compositional gradient index under the same lamination scheme and the span-to-thickness ratio leads to notable changes in the vertical displacements of Configuration 1 compared to Configuration 2. In Configuration 1, when examining 0°/0°, the vertical displacements exhibit 34.7% and 55.3% for {\lambda }_{\mathrm{s}}= 10, respectively. These values of α = 1 and 0.5 are 35.8% and 56.7% for {\lambda }_{\mathrm{s}}= 15; 36.1% and 57.2% for {\lambda }_{\mathrm{s}}= 20, respectively. Also, when observing 0°/90°, the vertical displacements show 26.3% and 46.1% decrease for {\lambda }_{\mathrm{s}}= 10, respectively. These values are 26.7% and 46.8% for {\lambda }_{\mathrm{s}}= 15; and 26.8% and 47.1% for {\lambda }_{\mathrm{s}}= 20, respectively. In Configuration 2, when assessing 0°/0°, the vertical displacements exhibit 17.3% and 27.8% for {\lambda }_{\mathrm{s}}= 10, respectively. These values of α = 1 and 0.5 are 16.8% and 26.6% for {\lambda }_{\mathrm{s}}= 15; 16.9% and 26.1% for {\lambda }_{\mathrm{s}}= 20, respectively. Also, when investigating 0°/90°, the vertical displacements show 16.0% and 26.8%{\lambda }_{\mathrm{s}}= 10, respectively. These values are 15.8% and 26.1% for {\lambda }_{\mathrm{s}}= 15; and 15.6% and 25.7% for {\lambda }_{\mathrm{s}}= 20, respectively. When variations in the lamination scheme are considered for each configuration and span-to-thickness ratio, the symmetric 0°/0° scheme exhibits the highest percentage differences, whereas the antisymmetric 0°/90° scheme shows the lowest. This behavior is attributed to the reduction in effective material properties in the 90° ply orientation of the x-axis of the beam, which decreases the changes in vertical displacement for both configurations in response to variation in the compositional gradient index.

The effect of the configurations: Tab.16 reveals that Configuration 1 exhibits the highest absolute vertical displacement, while Configuration 2 demonstrates the lowest, regardless of span-to-thickness ratios, compositional gradient index, and lamination schemes. The results of Configuration 2 are normalized with respect to Configuration 1. The percent differences between Configurations 1 and 2 rise as the compositional gradient index increases, indicating that the most significant reduction in vertical displacement occurs at α = 2 for both lamination schemes. At lower values of the compositional gradient index, the percentage differences for the 0°/0° and 0°/90° lamination schemes are relatively similar. However, as the compositional gradient index increases, the influence of the lamination scheme becomes more pronounced, resulting in a growing disparity between the percentage differences for the 0°/0° and 0°/90° lamination schemes. When comparing Configuration 1 to Configuration 2, keeping all other parameters such as compositional gradient index, and lamination schemes constant, the most significant reduction in vertical displacements is observed in the span-to-thickness ratio {\lambda }_{\mathrm{s}}= 20, while the smallest reduction is found in the {\lambda }_{\mathrm{s}}= 10 case (Fig.9).

4.3.2 Normal stress

The normal stress at the midspan of the beam is investigated to evaluate the effect of the compositional gradient index, different configurations, and lamination schemes in a three-phase composite beam with two-layer Configuration 1 and 2, with 0°/0° and 0°/90° layups. The span-to-thickness ratio is {\lambda }_{\mathrm{s}}= 10. The non-dimensional normal stress of the beam is obtained by {\tilde{\sigma }}_x(x,\mathit{\text{z}})=\left(w{h}^2/q{L}^2\right){\sigma }_x(x, \mathit{\text{z}}). First, variation of non-dimensional normal stresses ({\tilde{\sigma }}_x) at x={\displaystyle \frac{L}{2}} are plotted against the normalized height ratio and presented in Fig.10 for different compositional gradient indexes. For detailed analysis, the maximum tensile and compressive normal stresses ({\tilde{\sigma }}_x) through the thickness at the midspan of the beam for both Configuration 1 and Configuration 2 with symmetric and antisymmetric lamination schemes as follows.

The effect of the compositional gradient index (α): in Configuration 1 with a symmetric (0°/0°) lamination scheme, the evaluation of the normal stress distributions shows that the maximum tensile and compressive stress magnitudes remain nearly the same in the upper and lower regions of the cross-section (Fig.10(a)). However, as the composition gradient index α increases, the position of the absolute maximum stress values shows a slight shift toward the center of the cross-section. In contrast, for Configuration 1 with the antisymmetric (0°/90°) lamination scheme, the maximum tensile and compressive stress values are significantly influenced by changes in the compositional gradient index α, unlike the symmetric (0°/0°) lamination scheme (Fig.10(b)). The maximum tensile stress in the (0°/90°) configuration occurs in a region similar to that of the 0°/0° configuration. However, the maximum compressive stress reaches its extreme value in the upper part of the cross-section. This behavior is attributed to the reduction in effective material properties in the 90° ply orientation along the x-axis of the beam compared to the 0° ply orientation. In Configuration 2, whether employing a symmetric (0°/0°) or antisymmetric (0°/90°) lamination scheme (Fig.10(c) and Fig.10(d)), the maximum tensile and compressive stresses are located at the top and bottom regions of the cross-section, respectively. To analyze the effect of the compositional gradient index in detail, the maximum tensile and compressive normal stress values for α = 0.5 and α = 1 are normalized with respect to the corresponding results for α = 2 in each configuration and lamination scheme. In Configuration 1 with a symmetric (0°/0°) lamination scheme, the tensile or compressive normal stress of the α = 1 and 0.5 cases have –2.15% and –0.22%, respectively. Also, with antisymmetric (0°/90°) lamination scheme, the tensile normal stress of the α = 1 and 0.5 cases have –21.05% and –12.54%, respectively. The compressive normal stress of the α = 1 and 0.5 cases have 24.78% and 20.15%, respectively. In Configuration 2 with symmetric (0°/0°) lamination scheme, the tensile or compressive normal stress of the α = 1 and 0.5 cases have 19.12% and 24.70%, respectively. Also, with antisymmetric (0°/90°) lamination scheme, the tensile normal stress of the α = 1 and 0.5 cases have 17.46% and 22.99%, respectively. The compressive normal stress of the α = 1 and 0.5 cases have 17.10% and 25.08%, respectively.

The effect of the configurations: to analyze the effect of configuration in detail, the maximum tensile and compressive normal stresses of Configuration 1 are normalized with respect to those of Configuration 2 for each lamination scheme and compositional gradient index. For the symmetric (0°/0°) lamination scheme, the maximum percentage difference in tensile and compressive stress is 46.8%. For the antisymmetric (0°/90°) lamination scheme, the maximum percentage differences in tensile and compressive stresses are 55.1% and 47.3%, respectively. According to the percent differences, the antisymmetric (0°/90°) lamination scheme shows greater sensitivity to the change in ply orientation, as the lower stiffness of the 90° ply leads to significant variations in stress distribution. Tensile stress is particularly affected by this change, resulting in higher percent differences between Configuration 1 and Configuration 2. Conversely, the symmetric (0°/0°) lamination scheme exhibits smaller percent differences due to the higher stiffness of the 0° ply and its uniform orientation, resulting in a more consistent stress distribution.

4.3.3 Shear stress

The shear stress ({\tau }_{x\mathit{\text{z}}}) at the x=L/10 of the beam is investigated to assess the effect of the compositional gradient index, different configurations, and lamination schemes in a three-phase composite beam with a two-layer Configuration 1 and 2, with 0°/0° and 0°/90° layups. The span-to-thickness ratio is {\lambda }_{\mathrm{s}}= 10. The non-dimensional shear stress of the beam is obtained by {\tilde{\tau }}_{x\mathit{\text{z}}}(x ,\mathit{\text{z}})=(h w/ qL){\tau }_{x\mathit{\text{z}}} (x,\mathit{\text{z}}). First, variation of non-dimensional shear stresses \left({\tilde{\tau }}_{x\mathit{\text{z}}}\right) at x=L/10 are plotted against the normalized height ratio and presented in Fig.11 for different compositional gradient indexes. For detailed analysis, the absolute maximum shear stresses \left({\tau }_{x\mathit{\text{z}}}\right) through the thickness at the x=L/10 of the beam for both Configuration 1 and Configuration 2 with symmetric and antisymmetric lamination scheme as follows:

The effect of the compositional gradient index (\alpha): in both symmetric (0°/0°) and antisymmetric (0°/90°) lamination cases in Fig.11, as the compositional gradient index (α) increases, the absolute maximum shear stress values for Configuration 1 increase, while for Configuration 2, the absolute maximum shear stress values decrease. When the absolute maximum shear stress values for α = 0.5 and 1 are normalized with respect to the results for α = 2, for each configuration and lamination scheme. In Configuration 1, the shear stress shows 4.45% and 10.59% for symmetric (0°/0°) lamination scheme, respectively. These values of α = 1 and 0.5 are 0.09% and 3.57% for the antisymmetric (0°/90°) lamination scheme, respectively. In Configuration 2, the shear stress shows –4.41% and –8.69% for symmetric (0°/0°) lamination scheme, respectively. These values of α = 1 and 0.5 are –2.54% and –5.69% for the antisymmetric (0°/90°) lamination scheme, respectively. Based on the percentage differences, when examining the variation of the maximum shear stress values with respect to α, the change in the symmetric (0°/0°) lamination scheme is more pronounced compared to the change in the antisymmetric (0°/90°) lamination scheme.

The effect of the configurations: when the absolute maximum shear stresses of Configuration 1 are normalized with respect to Configuration 2. When comparing Configuration 1 to Configuration 2, keeping all other parameters constant, the most significant change in absolute maximum shear stress is observed in the compositional gradient index α = 2, while the smallest increase is found in the α = 0.5 case. When examining the percentage differences in absolute maximum shear stress values, it is observed that the percent differences of the symmetric (0°/0°) and antisymmetric (0°/90°) lamination scheme cases are relatively close to each other, with values as –19.1% and –21.0% at α = 0.5, respectively. However, for higher values of the α, the percent differences of the symmetric (0°/0°) and antisymmetric (0°/90°) lamination schemes increasingly diverge, with values as –44.8% and –32.7% at α = 2, respectively. It is concluded that the configuration effect becomes increasingly significant for higher values of the compositional gradient index, highlighting the influence of material property variations on stress distribution.

5 Conclusions

For the static and dynamic analysis of two-phase and multi-phase CNT-reinforced FG laminated composite beams, the precision performance of the presented mixed FE method inserted warping effect is about 0%–5.81% when compared with the three-dimensional solid elements. The W-MFE is based on Timoshenko beam theory and has 20-four DOFs. The principles of 3D elasticity theory are used to derive constitutive Eqs. (76, 95). W-MFE directly calculates forces and moments at nodal points, considers cross-sectional warping, offers high precision in stress calculations, time-saving computational advantage, and reduces the number of unknowns. Through the numerical examples presented in Subsubsection 4.1.1, the performance of the W-MFE results for the two-phase FG-CNTRC beam are thoroughly investigated and compared with the literature, including studies employing FSDT, SNDT, and higher order ZZT, as well as with three-dimensional solid element results (SOLID186). The key findings from the parametric analysis in Subsection 4.2 on free vibration and Subsection 4.3 on static analysis of three-phase CNT/polymer/fiber FG beams are summarized as follows: the free vibration analysis, due to the increase of compositional gradient index, the differences between Configuration 1, where the fiber-rich region is in the center of the cross-section, and Configuration 2, where the fiber-rich regions are located at the top and bottom of the cross-section, become more pronounced for each natural frequency. In static analysis, tensile and compressive stress distributions through the thickness reveal significant differences between configurations and lamination schemes. The peak tensile and compressive stress values in Configuration 1 are lower than those in Configuration 2, with these differences increasing as the compositional gradient index rises. Similarly, the peak shear stress values in Configuration 1 are lower than those in Configuration 2. However, increasing the compositional gradient index raises the peak shear stress in Configuration 1 while reducing it in Configuration 2. The compositional gradient index and fiber distribution within the cross-section can be tailored precisely to optimize the stress distribution, stiffness, displacement, and natural frequency characteristics of the three-phase CNT/polymer/fiber FG composite beams, depending on the desired structural requirements. The benchmark examples and results are presented in the literature for future reference and validation.

References

[1] Turaka S, Chintalapudi R, Kannaiyan Geetha N, Pappula B, Makgato S. Experimental and numerical analysis of the microstructure and mechanical properties of unidirectional glass fiber reinforced epoxy composites. Composite Structures, 2024, 331: 117887

[2] Taimour M, Abd-Elhady A A, Sallam H E D M, Sayed S A. Implementing functionally graded fibers technique to enhance pinned-joint performance in cross-ply laminate polymeric composites. Composite Structures, 2023, 313: 116931

[3] Pham Q H, Tran T T, Nguyen P C. Nonlinear dynamic analysis of functionally graded carbon nanotube-reinforced composite plates using MISQ20 element. Frontiers of Structural and Civil Engineering, 2023, 17(7): 1072–1085

[4] Altay U, Dorduncu M, Kadioglu S. An improved peridynamic approach for fatigue analysis of two dimensional functionally graded materials. Theoretical and Applied Fracture Mechanics, 2023, 128: 104152

[5] Dorduncu M, Apalak M K, Reddy J N. Stress wave propagation in a through-thickness functionally graded adhesive layer. Journal of Adhesion Science and Technology, 2019, 33(21): 2329–2355

[6] ErmisMOmurtag M H. Static and dynamic analysis of conical helices based on exact geometry via mixed FEM. International Journal of Mechanical Sciences, 2017, 131–132: 296–304

[7] SamaniegoEAnitescu CGoswamiSNguyen-ThanhV MGuoH HamdiaKZhuang XRabczukT. Rabczuk T. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112790

[8] GuoHZhuangX RabczukT. A deep collocation method for the bending analysis of Kirchhoff plate. 2021, arXiv: 2102.02617

[9] Zhuang X, Guo H, Alajlan N, Zhu H, Rabczuk T. Deep autoencoder based energy method for the bending, vibration, and buckling analysis of Kirchhoff plates with transfer learning. European Journal of Mechanics. A, Solids, 2021, 87: 104225

[10] Wang J, Duan X, Gong L, Nie S. Interfacial and filler size effects on mechanical/thermal/electrical properties of CNTs-reinforced nanocomposites. Polymers, 2024, 16(6): 808

[11] Ma H, Liu F, Wu Y, A X, Zhao Y. Axial compression tests and numerical simulation of steel reinforced recycled concrete short columns confined by carbon fiber reinforced plastics strips. Frontiers of Structural and Civil Engineering, 2022, 16(7): 817–842

[12] Feng X, Liu D, Guo Y, Zhong F, Zuo J, Liu W. Research on mechanical performance of longitudinal joints in segmental tunnel linings strengthened by fiber-reinforced plastic grid with polymer-cement-mortar method. Frontiers of Structural and Civil Engineering, 2024, 18(10): 1610–1625

[13] Irshidat M R, Al-Husban R S. Effect of bond enhancement using carbon nanotubes on flexural behavior of RC beams strengthened with externally bonded CFRP sheets. Frontiers of Structural and Civil Engineering, 2022, 16(1): 131–143

[14] Arda M, Aydogdu M. Vibration analysis of carbon nanotube mass sensors considering both inertia and stiffness of the detected mass. Mechanics Based Design of Structures and Machines, 2022, 50(3): 841–857

[15] Garg A, Chalak H D, Belarbi M O, Zenkour A M, Sahoo R. Estimation of carbon nanotubes and their applications as reinforcing composite materials—An engineering review. Composite Structures, 2021, 272: 114234

[16] Pektürk H Y, Demir B, Bilgi C, Öz F, Ersoy N. Effect of MWCNT on the properties of NCF-CFP hybrid composite fabricated via vacuum infusion. Journal of Reinforced Plastics and Composites, 2024, 43(21–22): 1214–1227

[17] Shokrieh M M, Rafiee R. On the tensile behavior of an embedded carbon nanotube in polymer matrix with non-bonded interphase region. Composite Structures, 2010, 92(3): 647–652

[18] Gok S G, Sengul O. Mechanical properties of alkali-activated slag based SIFCON incorporating waste steel fibers and waste glass. Construction and Building Materials, 2023, 408: 133697

[19] Özkılıç Y O, Beskopylny A N, Stel’makh S A, Shcherban E M, Mailyan L R, Meskhi B, Chernil’nik A, Ananova O, Aksoylu C, Madenci E. Lightweight expanded-clay fiber concrete with improved characteristics reinforced with short natural fibers. Case studies in construction materials, 2023, 19: e02367

[20] Khatir A, Capozucca R, Khatir S, Magagnini E, Benaissa B, Cuong-Le T. An efficient improved Gradient Boosting for strain prediction in Near-Surface Mounted fiber-reinforced polymer strengthened reinforced concrete beam. Frontiers of Structural and Civil Engineering, 2024, 18(8): 1148–1168

[21] Asyraf M R M, Hazrati K Z, Sheng D D C V, Rafidah M, Ashraf W, Madenci E, Özkılıç Y O, Aksoylu C, Nurazzi N M. Lignocellulosic sugar palm fibre-reinforced thermoplastic composites: Mechanical, thermal and dynamic mechanical properties. Fibers and Polymers, 2023, 24(8): 2625–2639

[22] Shokrieh M M, Rafiee R. Development of a full range multi-scale model to obtain elastic properties of CNT/polymer composites. Iranian Polymer Journal, 2012, 21(6): 397–402

[23] Rafiee R, Firouzbakht V. Multi-scale modeling of carbon nanotube reinforced polymers using irregular tessellation technique. Mechanics of Materials, 2014, 78: 74–84

[24] RafieeR. Carbon Nanotube-Reinforced Polymers: From Nanoscale to Macroscale. Amsterdam: Elsevier, 2017

[25] Rafiee R. Influence of carbon nanotube waviness on the stiffness reduction of CNT/polymer composites. Composite Structures, 2013, 97: 304–309

[26] Rafiee R, Zehtabzadeh H, Amini M R. Predicting mechanical properties of 3D printed nanocomposites using multi-scale modeling. Additive Manufacturing, 2024, 83(11): 104055

[27] Rafiee R, Zehtabzadeh H. Predicting the strength of carbon nanotube reinforced polymers using stochastic bottom-up modeling. Applied Physics. A, Materials Science and Processing, 2020, 126(8): 595

[28] RafieeRFirouzbakht V. Stochastic Multiscale Modeling of CNT/Polymer. Amsterdam: Elsevier, 2018, 503–520

[29] Rafiee R, Amohaji H. On the application of genetic algorithm for predicting the strength of CNT/ABS filaments using multi-scale modeling. Applied Physics. A, Materials Science and Processing, 2024, 130(7): 513

[30] RafieeRShokrieh M M. Prediction of mechanical properties of CNT based composites using multi-scale modeling and stochastic analysis. In: Tserpes K, Silvestre N, eds. Modeling of Carbon Nanotubes, Graphene and their Composites. Cham: Springer, 2014, 201–238

[31] Rafiee R, Shahzadi R. Predicting mechanical properties of nanoclay/polymer composites using stochastic approach. Composites. Part B, Engineering, 2018, 152: 31–42

[32] Rafiee R, Eskandariyun A. Predicting Young’s modulus of agglomerated graphene/polymer using multi-scale modeling. Composite Structures, 2020, 245: 112324

[33] Lan Y, Wang Y, Ren Z F. Physics and applications of aligned carbon nanotubes. Advances in Physics, 2011, 60(4): 553–678

[34] Marconnet A M, Yamamoto N, Panzer M A, Wardle B L, Goodson K E. Thermal conduction in aligned carbon nanotube–polymer nanocomposites with high packing density. ACS Nano, 2011, 5(6): 4818–4825

[35] Uzun B, Yaylı M Ö. Free vibration of a carbon nanotube-reinforced nanowire/nanobeam with movable ends. Journal of Vibration Engineering and Technologies, 2024, 12(4): 6847–6863

[36] Chowdhury S C, Haque B G, Okabe T, Gillespie J W Jr. Modeling the effect of statistical variations in length and diameter of randomly oriented CNTs on the properties of CNT reinforced nanocomposites. Composites. Part B, Engineering, 2012, 43(4): 1756–1762

[37] Mecklenburg M, Mizushima D, Ohtake N, Bauhofer W, Fiedler B, Schulte K. On the manufacturing and electrical and mechanical properties of ultra-high wt. % fraction aligned MWCNT and randomly oriented CNT epoxy composites. Carbon, 2015, 91: 275–290

[38] Subramanian P. Dynamic analysis of laminated composite beams using higher order theories and finite elements. Composite Structures, 2006, 73(3): 342–353

[39] Aribas U N. Bi-directional higher-order shear deformable mixed finite element formulation including couple effects for stresses of functionally graded curved 3D beams. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2024, 46(11): 641

[40] Bab Y, Kutlu A. Stress analysis of laminated HSDT beams considering bending extension coupling. Turkish Journal of Civil Engineering, 2023, 34(1): 1–23

[41] Sorrenti M, Gherlone M. Numerical and experimental predictions of the static behaviour of thick sandwich beams using a mixed {3,2}-RZT formulation. Finite Elements in Analysis and Design, 2024, 242: 104267

[42] Ermis M, Dorduncu M, Kutlu A. Peridynamic differential operator for stress analysis of imperfect functionally graded porous sandwich beams based on refined zigzag theory. Applied Mathematical Modelling, 2024, 133: 414–435

[43] Yurtsever B, Bab Y, Kutlu A, Dorduncu M. A new C0 continuous refined zigzag {1,2} finite element formulation for flexural and free vibration analyses of laminated composite beams. Composite Structures, 2024, 331: 117890

[44] Heshmati M, Yas M H. Free vibration analysis of functionally graded CNT-reinforced nanocomposite beam using Eshelby-Mori-Tanaka approach. Journal of Mechanical Science and Technology, 2013, 27(11): 3403–3408

[45] Yas M H, Heshmati M. Dynamic analysis of functionally graded nanocomposite beams reinforced by randomly oriented carbon nanotube under the action of moving load. Applied Mathematical Modelling, 2012, 36(4): 1371–1394

[46] Ansari R, Faghih Shojaei M, Mohammadi V, Gholami R, Sadeghi F. Nonlinear forced vibration analysis of functionally graded carbon nanotube-reinforced composite Timoshenko beams. Composite Structures, 2014, 113: 316–327

[47] Jam J E, Kiani Y. Low velocity impact response of functionally graded carbon nanotube reinforced composite beams in thermal environment. Composite Structures, 2015, 132: 35–43

[48] Kamarian S, Shakeri M, Yas M H, Bodaghi M, Pourasghar A. Free vibration analysis of functionally graded nanocomposite sandwich beams resting on Pasternak foundation by considering the agglomeration effect of CNTs. Journal of Sandwich Structures and Materials, 2015, 17(6): 632–665

[49] Ebrahimi F, Farazamandnia N. Thermo-mechanical analysis of carbon nanotube-reinforced composite sandwich beams. Coupled Systems Mechanics, 2017, 6: 207–227

[50] Borjalilou V, Taati E, Ahmadian M T. Bending, buckling and free vibration of nonlocal FG-carbon nanotube-reinforced composite nanobeams: exact solutions. SN Applied Sciences, 2019, 1(11): 1323

[51] Talebi S, Arvin H, Beni Y T. Thermal free vibration examination of sandwich piezoelectric agglomerated randomly oriented CNTRC Timoshenko beams regarding pyroelectricity. Engineering Analysis with Boundary Elements, 2023, 146: 500–516

[52] Lin F, Xiang Y. Vibration of carbon nanotube reinforced composite beams based on the first and third order beam theories. Applied Mathematical Modelling, 2014, 38(15–16): 3741–3754

[53] Jedari Salami S. Extended high order sandwich panel theory for bending analysis of sandwich beams with carbon nanotube reinforced face sheets. Physica E, Low-Dimensional Systems and Nanostructures, 2016, 76: 187–197

[54] Biswas S, Datta P. Finite element model for free vibration analyses of FG-CNT reinforced composite beams using refined shear deformation theories. IOP Conference Series. Materials Science and Engineering, 2021, 1206(1): 012019

[55] Madenci E. Free vibration analysis of carbon nanotube RC nanobeams with variational approaches. Advances in Nano Research, 2021, 11(2): 157–171

[56] Belarbi M O, Salami S J, Garg A, Daikh A A, Houari M S A, Dimitri R, Tornabene F. Mechanical behavior analysis of FG-CNT-reinforced polymer composite beams via a hyperbolic shear deformation theory. Continuum Mechanics and Thermodynamics, 2023, 35(2): 497–520

[57] Chalak H D, Zenkour A M, Garg A. Free vibration and modal stress analysis of FG-CNTRC beams under hygrothermal conditions using zigzag theory. Mechanics Based Design of Structures and Machines, 2023, 51(8): 4709–4730

[58] Nanoth R, Jayanarayanan K, Sarath Kumar P, Balachandran M, Pegoretti A. Static and dynamic mechanical properties of hybrid polymer composites: A comprehensive review of experimental, micromechanical and simulation approaches. Composites. Part A, Applied Science and Manufacturing, 2023, 174: 107741

[59] Madenci E, Özkılıç Y O, Aksoylu C, Asyraf M R M, Syamsir A, Supian A B M, Mamaev N. Buckling analysis of CNT-reinforced polymer composite beam using experimental and analytical methods. Materials, 2023, 16(2): 614

[60] Rafiee M, He X Q, Mareishi S, Liew K M. Modeling and stress analysis of smart CNTs/fiber/polymer multiscale composite plates. International Journal of Applied Mechanics, 2014, 6(3): 1450025

[61] Bacciocchi M, Tarantino A M. Natural frequency analysis of functionally graded orthotropic cross-ply plates based on the finite element method. Mathematical and Computational Applications, 2019, 24(2): 52

[62] Bacciocchi M, Tarantino A M. Time-dependent behavior of viscoelastic three-phase composite plates reinforced by Carbon nanotubes. Composite Structures, 2019, 216: 20–31

[63] Tornabene F, Bacciocchi M, Fantuzzi N, Reddy J N. Multiscale approach for three phase CNT/polymer/fiber laminated nanocomposite structures. Polymer Composites, 2019, 40(S1): E102–E126

[64] Bacciocchi M. Buckling analysis of three-phase CNT/polymer/fiber functionally graded orthotropic plates: Influence of the non-uniform distribution of the oriented fibers on the critical load. Engineering Structures, 2020, 223: 111176

[65] Zhu P, Xu J, Qu W. Fatigue shear performance of concrete beams reinforced with hybrid (glass-fiber-reinforced polymer + steel) rebars and stirrups. Frontiers of Structural and Civil Engineering, 2021, 15(3): 576–594

[66] Rahmani S, Seidi H, Kakavand F, Soltani E. Free vibration analysis of a cantilever trapezoidal sandwich plate with an auxetic reentrant honeycomb core and laminated nanocomposite polymer/CNT/fiber face sheets. Mechanics Based Design of Structures and Machines, 2024, 52(9): 7180–7210

[67] ZhangSZhang YLeiQYangYWangY XuFYanZ ZhuH. Zhu H. Acoustic emission characteristics of damage evolution of multi-scale fiber reinforced rubberized concrete under uniaxial compression and tension after being subjected to high temperatures. Frontiers of Structural and Civil Engineering, 2024, 18(8): 1237–1266

[68] He X Q, Rafiee M, Mareishi S, Liew K M. Large amplitude vibration of fractionally damped viscoelastic CNTs/fiber/polymer multiscale composite beams. Composite Structures, 2015, 131: 1111–1123

[69] Afshin S, Yas M H. Dynamic and buckling analysis of polymer hybrid composite beam with variable thickness. Applied Mathematics and Mechanics, 2020, 41(5): 785–804

[70] Alambeigi K, Mohammadimehr M, Bamdad M, Rabczuk T. Free and forced vibration analysis of a sandwich beam considering porous core and SMA hybrid composite face layers on Vlasov’s foundation. Acta Mechanica, 2020, 231(8): 3199–3218

[71] Dabbagh A, Rastgoo A, Ebrahimi F. Thermal buckling analysis of agglomerated multiscale hybrid nanocomposites via a refined beam theory. Mechanics Based Design of Structures and Machines, 2021, 49(3): 403–429

[72] Daikh A A, Belarbi M O, Salami S J, Ladmek M, Belkacem A, Houari M S A, Ahmed H M, Eltaher M A. Eltaher M A. A three-unknown refined shear beam model for the bending of randomly oriented FG-CNT/fiber-reinforced composite laminated beams rested on a new variable elastic foundation. Acta Mechanica, 2023, 234(10): 5171–5186

[73] Li R, Ding R, Hozuri A. Nonlinear vibrations of multiscale composite beams on a nonlinear softening foundation. European Journal of Mechanics. A, Solids, 2023, 100: 104993

[74] Aribas U N, Ermis M, Omurtag M H. The static and stress analyses of axially functionally graded exact super-elliptical beams via mixed FEM. Archive of Applied Mechanics, 2021, 91(12): 4783–4796

[75] Aribas U N, Aydin M, Atalay M, Omurtag M H. Cross-sectional warping and precision of the first-order shear deformation theory for vibrations of transversely functionally graded curved beams. Applied Mathematics and Mechanics, 2023, 44(12): 2109–2138

[76] Ermis M. Warping-included mixed FE approach of beating characteristics in functionally graded graphene platelet-reinforced composite spatially curved beams under harmonic excitation force. Archive of Applied Mechanics, 2024, 94(12): 3687–3713

[77] Yas M H, Samadi N. Free vibrations and buckling analysis of carbon nanotube-reinforced composite Timoshenko beams on elastic foundation. International Journal of Pressure Vessels and Piping, 2012, 98: 119–128

[78] Vo-Duy T, Ho-Huu V, Nguyen-Thoi T. Free vibration analysis of laminated FG-CNT reinforced composite beams using finite element method. Frontiers of Structural and Civil Engineering, 2019, 13(2): 324–336

[79] Karamanli A, Vo T P. Finite element model for carbon nanotube-reinforced and graphene nanoplatelet-reinforced composite beams. Composite Structures, 2021, 264: 113739

[80] Garg A, Chalak H D, Zenkour A M, Belarbi M O, Sahoo R. Bending and free vibration analysis of symmetric and unsymmetric functionally graded CNT reinforced sandwich beams containing softcore. Thin-walled Structures, 2022, 170: 108626

[81] Kumar P, Srinivas J. Free vibration, bending and buckling of a FG-CNT reinforced composite beam: Comparative analysis with hybrid laminated composite beam. Multidiscipline Modeling in Materials and Structures, 2017, 13(4): 590–611

[82] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica, 1973, 21(5): 571–574

[83] Shi X, Hassanzadeh-Aghdam M K, Ansari R. A comprehensive micromechanical analysis of the thermoelastic properties of polymer nanocomposites containing carbon nanotubes with fully random microstructures. Mechanics of Advanced Materials and Structures, 2019, 28(4): 331–342

[84] Tornabene F, Fantuzzi N, Bacciocchi M, Viola E. Effect of agglomeration on the natural frequencies of functionally graded carbon nanotube-reinforced laminated composite doubly-curved shells. Composites. Part B, Engineering, 2016, 89: 187–218

[85] Dabbagh A, Rastgoo A, Ebrahimi F. Post-buckling analysis of imperfect multi-scale hybrid nanocomposite beams rested on a nonlinear stiff substrate. Engineering with Computers, 2022, 38(1): 301–314

[86] Hahn H T. Simplified formulas for elastic moduli of unidirectional continuous fiber composites. Journal of Composites Technology and Research, 1980, 2(3): 5–7

[87] Shi D L, Feng X Q, Huang Y Y, Hwang K C, Gao H. The effect of nanotube waviness and agglomeration on the elastic property of carbon nanotube-reinforced composites. Journal of Engineering Materials and Technology, 2004, 126(3): 250–257

[88] JonesR M. Mechanics of Composite Materials. Boca Raton: CRC Press, 1999

[89] Yıldırım V. Governing equations of initially twisted elastic space rods made of laminated composite materials. International Journal of Engineering Science, 1999, 37(8): 1007–1035

[90] Aribas U N, Ermis M, Eratli N, Omurtag M H. The static and dynamic analyses of warping included composite exact conical helix by mixed FEM. Composites. Part B, Engineering, 2019, 160: 285–297

[91] Yousefi A, Rastgoo A. Free vibration of functionally graded spatial curved beams. Composite Structures, 2011, 93(11): 3048–3056

[92] Bhimaraddi A, Chandrashekhara K. Some observations on the modeling of laminated composite beams with general layups. Composite Structures, 1991, 19(4): 371–380

[93] Aribas U N, Ermis M, Kutlu A, Eratli N, Omurtag M H. Forced vibration analysis of composite-geometrically exact elliptical cone helices via mixed FEM. Mechanics of Advanced Materials and Structures, 2022, 29(10): 1456–1474

[94] Jog C S, Mokashi I S. A finite element method for the Saint-Venant torsion and bending problems for prismatic beams. Computers and Structures, 2014, 135: 62–72

[95] Aribas U N, Atalay M, Omurtag M H. Warping included mixed finite elements for bending and stresses of functionally graded exact curved beams. Mechanics of Advanced Materials and Structures, 2024, 31(20): 5040–5056

[96] Omurtag M H, Aköz A Y. The mixed finite element solution of helical beams with variable cross-section under arbitrary loading. Computers and Structures, 1992, 43(2): 325–331

[97] OdenJ TReddy J N. Variational Methods in Theoretical Mechanics. Berlin: Springer, 2012

[98] Doğruoğlu A N, Omurtag M H. Stability analysis of composite-plate foundation interaction by mixed FEM. Journal of Engineering Mechanics, 2000, 126(9): 928–936

[99] Omurtag M H, Aköz A Y. A compatible cylindrical shell element for stiffened cylindrical shells in a mixed finite element formulation. Computers and Structures, 1993, 49(2): 363–370

[100] Omurtag M H, Akoz A Y. Hyperbolic paraboloid shell analysis via mixed finite element formulation. International Journal for Numerical Methods in Engineering, 1994, 37(18): 3037–3056

[101] Eratli N, Yilmaz M, Darilmaz K, Omurtag M H. Dynamic analysis of helicoidal bars with non-circular cross-sections via mixed FEM. Structural Engineering and Mechanics, 2016, 57(2): 221–238

[102] ANSYS®Academic Research Mechanical. Release 17.1. Canonsburg, PA: ANSYS, Inc. 2017

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn