*Article* **A Multi-Scale Method for Designing Hybrid Fiber-Reinforced Composite Drive Shafts with Carbon Nanotube Inclusions**

**Stelios K. Georgantzinos 1,2,\*, Panagiotis A. Antoniou <sup>2</sup> and Stylianos I. Markolefas <sup>2</sup>**


**Abstract:** In this paper, the modal and linear buckling analysis of a laminated composite drive shaft reinforced by 11 multi-walled carbon nanotubes (MWCNTs) was carried out using an analytical approach, as well as the finite element method (FEM). The theoretical model is based on classical laminated theory (CLT). The fundamental frequency and the critical buckling torque were determined for different fiber orientation angles. The Halpin–Tsai model was employed to calculate the elastic modulus of composites having randomly oriented nanotubes. The effect of various carbon nanotube (CNT) volume fractions in the epoxy resin matrix on the material properties of unidirectional composite laminas was also analyzed. The fundamental frequency and the critical buckling torque obtained by the finite element analysis and the analytical method for different fiber orientation angles were in good agreement with each other. The results were verified with data available in the open literature, where possible. For the first time in the literature, the influence of CNT fillers on various composite drive shaft design parameters such as the fundamental frequency, critical speed, and critical buckling torque of a hybrid fiber-reinforced composite drive shaft is finally predicted.

**Keywords:** buckling; fundamental frequency; finite element method; laminated composites; drive shaft; carbon nanotubes

### **1. Introduction**

Drive shafts play important roles in the transmission vehicle system. They are used to transmit motion from the differential to the wheels. As a result, torsional, bending, and normal forces occur during its operation [1,2]. The demand for more powerful engines with higher torques and the new domains of application with increased stresses, as well as the improved materials and the new production processes, have imposed completely new operational requirements on drive shafts as a mechanical component [3–5]. Meanwhile, major progress has been made from the automotive industry in the field of drive shafts design, using new manufacturing processes, to conform to the constantly growing demands of various customers [6–11].

Over recent years, a rapid development has been observed concerning drive shafts made of composite materials that, in addition to their low weight, are able to respond effectively to their functional demands. Shinde et al. [12] designed a glass-epoxy composite drive shaft for light motor vehicles and investigated it for torsional strength, natural frequency, and critical speed. The results obtained were compared with a steel drive shaft. The developed composite shaft has proved to be the best alternative to the conventional steel shaft as its replacement resulted in a 73% weight reduction of the drive shaft. Nadeem et al. [13] carried out a review on the design and analysis of the performance of composite drive shafts made of different materials, such as carbon, glass, Kevlar, and boron with epoxy resin. They concluded that the fiber orientation angle and stacking sequence have great influences on the buckling torque and on the dynamic characteristics. Yefa et al. [14]

**Citation:** Georgantzinos, S.K.; Antoniou, P.A.; Markolefas, S.I. A Multi-Scale Method for Designing Hybrid Fiber-Reinforced Composite Drive Shafts with Carbon Nanotube Inclusions. *J. Compos. Sci.* **2021**, *5*, 157. https://doi.org/10.3390/jcs5060157

Academic Editor: Francesco Tornabene

Received: 29 May 2021 Accepted: 8 June 2021 Published: 10 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

implemented a new mechanical analytical solution of torsional stiffness for a composite drive shaft with balance laminate, based on CLT and mechanical analysis. Finite element analysis, CLT, and experiments were carried for carbon-fiber-reinforced polymers (CFRP) drive shafts to predict the torsional stiffness. Other valuable works on both the FEM and CFRP composite materials can be found in [15,16].

Fiber-reinforced polymer matrix composites have emerged as a major class of structural materials and are being considered for use as a substitute for many traditional metallic materials in a large number of weight-critical components in aerospace, automotive, and other industries, due to their high strength/weight ratio and stiffness/weight ratio [17,18]. The simultaneous development of nanotechnology has given a new perspective on the use of composite materials, as the enhanced characteristics of nanostructures are able to improve significantly the already important mechanical properties of classic composite materials.

Nanocomposites have the potential of becoming the future structural material not only in terms of greater mechanical properties, but also in terms of superior thermal, electrical, optical, and other properties [19–22]. Polymer nanocomposites are polymer matrix composites in which the reinforcement has at least one of its dimensions in the nanometer range. The mechanical properties of the nano-reinforcement are considerably high and the ratio of their surface area to volume is high as well, which means that a great interfacial interaction with the matrix can be provided. These are the leading reasons for such highly improved properties in nanocomposites [23]. The usage of polymer nanocomposites as a material for primary structures is in its early stage, but their potential in future aircraft structures has been realized [24]. Lockheed Martin had revealed that the F-35 Lightning II will be the first mass-produced aircraft to integrate structural nanocomposites in nonload-bearing airframe components. A thermoset epoxy reinforced by CNTs will replace carbon fiber as the material used to produce F-35 wingtip fairings, beginning with low-rate initial production (LRIP)-4 aircraft [25]. Shakil et al. [26] carried out a review on the properties and fabrication techniques of fiber-reinforced polymer nanocomposites subjected to simulated accidental ballistic impacts. High-specific-surface-area nanoparticles were used for matrix modification to induce nano-scale toughness mechanisms, with a focus on the ballistic performance of composite structures. Tüzemen et al. [27] investigated the effects of nanoclay (NC), CNT, and a hybrid of both inclusions on the bending, tensile, and bearing strengths of nanocomposite plates. Mechanical properties were investigated by applying bending and tensile testing, as well as tensile testing, with bolted joints on the nanocomposite plates.

Khoramishad et al. [28] investigated the effect of adding MWCNTs on the high-velocity impact behavior of fiber metal laminates (FMLs). They concluded that incorporating 0.5 wt% of MWCNTs into the composite laminate of FML resulted in the maximum reduction of 29.8% in projectile residual velocity and the maximum increase of 18.9% in the absorbed energy during projectile perforation compared to the unreinforced FMLs. Chiou et al. [29] studied the synergistic effects on the mechanical properties of composites materials by adding graphene nanoplatelets (GNPs) and nanocarbon aerogel (NCA) in the epoxy resins. The experimental results demonstrated that the mechanical properties of GNP/NCA/epoxy nanocomposites and GNP/NCA/CFRP laminates were optimized via reinforcement upon the addition of GNP/NCA hybrids.

In the last two decades, several exotic forms of low-dimensional carbon inclusions have been discovered, including fullerenes (zero dimension), CNTs (one dimension), and graphene (two dimension). All these materials share graphite's sp2 carbon bonding. Graphene, the mother of all graphitic forms, is a monolayer of carbon atoms tightly packed into a two-dimensional hexagonal lattice sheet that can be wrapped up to fullerenes, rolled into CNTs, and stacked into graphite [30]. Since the discovery of CNTs by Iijima in 1991 [31], research on their growth, characterization, and application development has exploded [32,33]. Their extraordinary physical properties have drawn the interest of the scientific community studying for numerous applications in nanoelectromechanical systems and nanoelectronics [34].

Ahmadipour et al. [35] studied the synthesis of calcium copper titanate (CaCu3Ti4O12)/ MWCNT composites, using the ultrasonic technique. The effects of MWCNT content on the structural, dielectric, and mechanical properties of copper calcium titanate (CCTO) were investigated by dielectric measurement, as well as tensile strength and flexural strength tests. An important observation was that the CCTO mixed with 0.2 wt% MWCNT exhibited the highest dielectric permittivity and the lowest dielectric loss. With the addition of 0.2 wt% MWCNT, the values of load, tensile, and flexural strength increased to 10.38 kN, 101.88 MPa, and 275.07 MPa, respectively, due to the improvement in densification. These results have values for the fabrication of CCTO and the optimization of its performance for electronic devices such as capacitors and antennas.

Recently, a lot of attention has been also given by researchers to graphene nanosheets (GNS) and CNT/polymer matrix composite materials. Georgantzinos et al. studied the mechanical elastic behavior [36,37] and vibration response [38] of laminated polymer composite plates with carbon nanostructure inclusions, using a multilevel framework, starting from the nanoscale, up to the laminated hybrid composite plates. Rafiee et al. [39] compared the mechanical properties of epoxy resin nanocomposites with GNS, single-walled CNTs, and multi-walled CNTs at the same nanofiller weight fraction. In an excellent recent work, Tas and Soykok [40] determined theoretically the engineering constants of CNTbased composite laminas. Bending analysis was performed on a composite plate under concentrated and distributed load. The results showed that elastic constants increased with the added CNT fraction. Additionally, the flexural modulus of the laminated composite plate showed significant improvement and the maximum deflection decreased.

Various techniques and approaches on the modeling of the mechanical properties of composites and nanostructures have been developed. Stamoulis et al. [41] dealt with the numerical modeling of the low-velocity impact damage of laminated composites, which have increasingly important applications in aerospace primary structures. Their purpose was to present and validate a computationally efficient approach in order to explore the effect of critical parameters on the impact damage characteristics. Barretta et al. [42] focused on modeling the bending of armchair CNTs by means of gradient elasticity theory. The estimation of Young's modulus depending on the armchair CNT diameter and length was investigated. According to the results presented, Young's modulus was found to be dependent on small size effects and ranged between 816.87 and 1067.24 MPa.

Georgantzinos et al. [43] investigated the thermomechanical buckling of single-walled CNTs (SWCNTs) by a structural mechanics method, providing theoretical predictions, concerning the compressive buckling response of SWCNTs of different chiralities and sizes under thermal conditions. Aghdam et al. [44] proposed a new micromechanics approach for predicting the elastic modulus of randomly oriented and distributed wavy CNTsreinforced polymer nanocomposites. The results revealed that the proposed method has a good accuracy according to the experimental data as compared. Mehar and Panda [45] investigated the deflection behavior of a composite plate reinforced with CNTs, numerically using FEM, and the result accuracy was verified via three-point experimental bending test data. The Young's modulus for plates consisting of randomly oriented MWCNTs obtained by using the Mori–Tanaka method showed an increase with the volume fraction of MWCNTs. Georgantzinos [46,47] developed a novel computational model for the prediction of the mechanical behavior of graphene and graphyne monolayers, according to their atomistic structure. The harmonic approximation was utilized for describing the interaction potential energies based on molecular theory. Spanos et al. [48] described a micromechanical finite element approach for the estimation of the elastic mechanical properties of graphene-reinforced composites. The load transfer conditions between the graphene and the matrix were modeled using joint elements connecting the two materials, simulating the interfacial region.

Feng et al. [49] investigated the mechanical properties of a graphene/PMMA nanocomposite system by using the molecular dynamics simulations. The results showed that the Young's and shear moduli increased with the graphene volume fraction and decreased as

the temperature rose from 300 to 500 K, while the efficiency of the reinforcement reduced with the graphene content. Rui et al. [50] presented an atomistic study on the mechanical behaviors of a polymer nanocomposite reinforced with defective graphene using molecular dynamics simulations, with a particular focus on the influences of temperature changes and atom vacancy defects. The Halpin–Tsai model was then modified based on the results to enable the effects of temperature and graphene defects to be considered in determining the Young's modulus of graphene-reinforced nanocomposites.

Even though a comprehensive investigation, in general, concerning polymer nanocomposites has been performed in the open literature, as well as their possible applications in automotive industry already being noticed [51], the potential performance improvement of specific components due to the use of nanomaterials has not been extensively reported yet. In this study, the development of suitable computational procedures, based on theoretical and finite element analysis, for the prediction of the mechanical behavior of automotive drive shafts manufactured from laminated composite materials, reinforced with CNTs, was presented. To the author's best knowledge, this study is the first one in the open literature concerning the effect of CNTs on the mechanical performance of composite drive shafts, providing novel results. It aims to offer a relatively simple and complete procedure for the prediction of the critical design parameters helping the design process and development of drive shafts of enhanced performance.

### **2. Theoretical Approach**

### *2.1. Unidirectional Composite Lamina with CNT Inclusions*

The first step to determine the mechanical response of nanocomposites under certain loads is the estimation of their material properties such as Young's modulus, Poisson's ratio, and the shear modulus. These elastic constants of nano-reinforced composite laminas may be calculated theoretically using the combination of simple models such as the Halpin– Tsai model, the Tsai–Pagano equation, and the Rule of mixtures. The models can be implemented considering certain assumptions. Here, CNTs were considered as randomly and homogeneously oriented into the matrix, while the CNTs-polymer (hybrid) matrix exhibited isotropic behavior. The unidirectional fiber-reinforced composite lamina with CNT inclusions is shown in Figure 1.

**Figure 1.** Schematic illustration of a unidirectional fiber-reinforced composite lamina, with the inclusion of CNTs in the polymer matrix. CNTs are considered as randomly oriented.

### *2.2. Hybrid CNT-Polymer Matrix Elastic Constants*

The well-established Halpin–Tsai model for fiber-reinforced composites was utilized to predict the Young's modulus of the nanocomposites. The CNTs were considered as randomly oriented long discontinuous fiber laminas. The Young's modulus of a CNTadded matrix can be calculated from the following equation [39]:

$$E\_{m-cut} = E\_m \left( \frac{3}{8} \frac{1 + 2\left(\frac{l\_{cut}}{d\_{cnt}}\right) \beta\_1 V\_{cnt}}{1 - \beta\_1 V\_{cnt}} + \frac{5}{8} \frac{1 + 2\beta\_2 V\_{cnt}}{1 - \beta\_2 V\_{cnt}} \right) \tag{1}$$

where

$$\beta\_1 = \frac{\left(\frac{E\_{eq}}{E\_m}\right) - 1}{\left(\frac{E\_{eq}}{E\_m}\right) + 2\left(\frac{I\_{cat}}{d\_{cat}}\right)} \text{ and } \beta\_2 = \frac{\left(\frac{E\_{eq}}{E\_m}\right) - 1}{\left(\frac{E\_{eq}}{E\_m}\right) + 2} \tag{2}$$

*Em* is the Young's modulus of the matrix, *Vcnt* is the volume content of CNTs, while *dcnt* and *lcnt* are the average outer diameter and length of the nanotube, respectively. *Eeq* is the equivalent modulus of nanotubes considering the hollow tube as a solid cylinder and can be expressed as:

$$E\_{cq} = \left(\frac{2t}{r}\right) E\_{cut} \tag{3}$$

where *t*, *r*, and *Ecnt* are the wall thickness, radius, and Young's modulus of the nanotube, respectively.

The shear modulus of the CNT-polymer matrix can be calculated by using the Tsai– Pagano equation [40]:

$$G\_{m-cnt} = E\_m \left( \frac{1}{8} \frac{1 + 2\left(\frac{l\_{cnt}}{d\_{cnt}}\right) \beta\_1 V\_{cnt}}{1 - \beta\_1 V\_{cnt}} + \frac{1}{4} \frac{1 + 2\beta\_2 V\_{cnt}}{1 - \beta\_2 V\_{cnt}} \right) \tag{4}$$

The Poisson's ratio for the CNT-polymer matrix exhibiting isotropic behavior can be calculated from the following equation:

$$v\_{m-cnt} = \left(\frac{E\_{m-cnt}}{2G\_{m-cnt}}\right) - 1\tag{5}$$

The density of the nanocomposite (hybrid) matrix can also be calculated from the following equation:

$$
\rho\_{m-\text{cut}} = \rho\_{\text{cut}} V\_{\text{cnt}} + \rho\_m V\_m \tag{6}
$$

where *ρcnt* is the density of the nanotube, *Vm* is the volume fraction of the matrix, and *ρ<sup>m</sup>* is the matrix density.

### *2.3. Unidirectional Composite Lamina Elastic Constants*

The material properties of the unidirectional composite lamina, adopting the rule of mixtures, may be calculated from the following equations [23,52]:

Longitudinal Young's modulus *E*<sup>1</sup> = *Ef Vf* + *Em*−*cntVm*−*cnt* (7)

$$\text{Transverse Young's modulus } E\_2 = \frac{E\_f E\_{m-cnt}}{E\_f V\_{m-cnt} + E\_{m-cnt} V\_f} \tag{8}$$

$$\text{In-plane Shear modulus } G\_{12} = \frac{G\_f G\_{m-cnt}}{G\_f V\_{m-cnt} + G\_{m-cnt} V\_f} \tag{9}$$

$$\text{Out-of-plane Shear modulus } G\_{23} = \frac{E\_2}{2(1+v\_{23})}\tag{10}$$

$$\text{Major Poisson's ratio } \upsilon\_{12} = \upsilon\_f V\_f + \upsilon\_{m-cnt} V\_{m-cnt} \tag{11}$$

$$\text{Transverse Poisson's ratio } v\_{23} = v\_{12} \frac{(1 - v\_{21})}{(1 - v\_{12})} \tag{12}$$

where *Ef*, *Vf*, *vf*, *Gf*, and *Vm-cnt* are Young's modulus of the fiber, the volume fraction of the fiber, Poisson's ratio of the fiber, the shear modulus of the fiber, and the volume fraction of the hybrid matrix, respectively.

Alternatively, semi-empirical models may be adopted for the evaluation of the material properties of the unidirectional composite lamina. An efficient model is the Halpin–Tsai [53] model, which can be used over a wide range of elastic properties and fiber volume fractions. In this approach, the longitudinal Young's modulus *E*<sup>1</sup> and major Poisson's ratio *ν*<sup>12</sup> are evaluated by Equations (7) and (11), respectively, while the other engineering constants are provided by the following equation:

$$\frac{P}{P\_{m-cnt}} = \left(\frac{1 + \xi \eta \, V\_f}{1 - \eta \, V\_f}\right) \tag{13}$$

where *Pm-cnt* are the related properties of the CNT-polymer matrix, and *P* can be considered as the transverse Young's modulus *E*2, transverse Poisson's ratio *ν*23, in-plane shear modulus *G*12, and out-of-plane shear modulus *G*23. *η* is an experimental factor computed by using the next expression:

$$\eta = \frac{\left(\frac{P\_f}{P\_{m-cnt}}\right) - 1}{\left(\frac{P\_f}{P\_{m-cnt}}\right) + \xi} \tag{14}$$

The value of the reinforcing factor *ξ* depends on the fiber geometry, packing geometry, and loading conditions. For circular fibers in a square array, *ξ* = 2 for *E*2, and *ξ* = 1 for *ν*23, *G*12, and *G*23, as referred in [52].

### *2.4. Theoretical Modeling of Composite Drive Shaft*

CLT was used to analyze the behavior of the composite laminated drive shaft. The following assumptions must be considered in order to make the problem solvable [54,55]:


Given the previous assumptions, a lamina can be analyzed as a plane stress problem and the stress–strain relations can take the following matrix form [55]:

$$\left\{ \begin{array}{c} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \gamma\_{12} \end{array} \right\} = [S] \left\{ \begin{array}{c} \sigma\_{11} \\ \sigma\_{22} \\ \tau\_{12} \end{array} \right\} \tag{15}$$

The square matrix on the right-hand side of Equation (15), denoted by [*S*], is called the reduced compliance matrix. The elements of [*S*] are given by:

$$\mathbf{S}\_{11} = \frac{1}{E\_1}, \mathbf{S}\_{12} = -\frac{\nu\_{12}}{E\_1}, \mathbf{S}\_{22} = \frac{1}{E\_2}, \mathbf{S}\_{66} = \frac{1}{G\_{12}}, \mathbf{S}\_{16} = \mathbf{S}\_{26} = 0\tag{16}$$

The inverse form of Equation (15) is as follows:

$$\left\{ \begin{array}{c} \sigma\_{11} \\ \sigma\_{22} \\ \tau\_{12} \end{array} \right\} = [Q] \left\{ \begin{array}{c} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \gamma\_{12} \end{array} \right\} \tag{17}$$

The square matrix on the right-hand side of Equation (17), denoted by [*Q*], is called the reduced stiffness matrix, in which:

$$[Q] = \begin{bmatrix} Q\_{11} & Q\_{16} & 0 \\ Q\_{12} & Q\_{22} & 0 \\ 0 & 0 & Q\_{66} \end{bmatrix} \tag{18}$$

and the elements of [*Q*] are given by:

$$Q\_{11} = \frac{E\_1}{1 - \nu\_{12}\nu\_{21}}, Q\_{12} = \frac{\nu\_{12}E\_2}{1 - \nu\_{12}\nu\_{21}}, Q\_{22} = \frac{E\_2}{1 - \nu\_{12}\nu\_{21}}, Q\_{66} = G\_{12}, \nu\_{21} = \frac{E\_2}{E\_1}\nu\_{12} \tag{19}$$

Unidirectional laminae are weak in transverse mechanical properties; to overcome this, they are usually placed in different orientations in a laminate. It is necessary to know the stress–strain relations for a generally orthotropic lamina [55].

We can relate global stresses to global strains as:

$$\left\{ \begin{array}{c} \sigma\_{11} \\ \sigma\_{22} \\ \tau\_{12} \end{array} \right\} = \left[T\right]^{-1} \left[Q\right] \left[R\right] \left[T\right] \left[R\right]^{-1} \left\{ \begin{array}{c} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \gamma\_{12} \end{array} \right\} \tag{20}$$

The matrix *Q* = [*T*] −1 [*Q*][*R*][*T*][*R*] <sup>−</sup><sup>1</sup> is called the transformed reduced stiffness matrix and it is obtained by multiplying the five 3 × 3 matrices on the right-hand side of Equation (20).

By carrying out the necessary multiplication and simplification, it can be shown that the elements of *Q* are given by:

$$\overline{Q\_{11}} = Q\_{11}\cos^4\theta + 2(Q\_{12} + 2Q\_{66})\sin^2\theta\cos^2\theta + Q\_{22}\sin^4\theta\tag{21}$$

$$\overline{Q\_{12}} = (Q\_{11} + Q\_{22} - 4Q\_{66})\sin^2\theta\cos^2\theta + Q\_{12}(\sin^4\theta\cos^4\theta) \tag{22}$$

$$\overline{Q\_{16}} = (Q\_{11} - Q\_{12} - 2Q\_{66})\sin\theta\cos^3\theta + (Q\_{12} - Q\_{22} + 2Q\_{66})\sin^3\theta\cos\theta\tag{23}$$

$$\overline{Q\_{22}} = Q\_{11}\sin^4\theta + 2(Q\_{12} + 2Q\_{\ell\theta})\sin^2\theta\cos^2\theta + Q\_{22}\cos^4\theta\tag{24}$$

$$\overline{Q\_{26}} = (Q\_{11} - Q\_{12} - 2Q\_{66})\sin^3\theta\cos\theta + \left(Q\_{12} - Q\_{22} + 2Q\_{66}\right)\sin\theta\cos^3\theta \tag{25}$$

$$\overline{Q\_{\theta\theta}} = (Q\_{11} + Q\_{22} - 2Q\_{12} - 2Q\_{\theta\theta})\sin^2\theta\cos^2\theta + Q\_{\theta\theta}(\sin^4\theta\cos^4\theta) \tag{26}$$

The elements of *Q* are functions of stiffness elements *Q*11, *Q*12, *Q*22, and *Q*<sup>66</sup> and the lamina angle *θ*. As the elements of [*Q*] are functions of four engineering constants, the elements of *Q* are consequently functions of four engineering constants and the lamina angle *θ* [55].

The next step is to construct the extensional stiffness matrix [*A*], which is the summation of the products of the transformed reduced stiffness matrix *Q* of the layer and the respective thicknesses [56].

$$A\_{ij} = \sum\_{k=1}^{n} \left\{ Q\_{ij} \right\}\_n (z\_k - z\_{k-1}) \tag{27}$$

where *i,j* = 1,2,6, the summation is carried over all *n* plies of the laminate, and *zk*, *zk*−<sup>1</sup> are the upper and lower z coordinates of the *kth* ply, respectively.

Equation (21) describes the membrane deformations of a laminate under in-plane loads and incorporates the Kirchhoff plate theory assumption that plane sections remain planar and perpendicular to the neutral axis.

The matrix [*A*] is used to calculate *Ex* and *Eh*, which are the average moduli in the axial and hoop directions, respectively, and *t* is the total thickness of drive shaft.

$$E\_{\rm x} = \frac{1}{t} \left[ A\_{11} - \frac{A\_{12}^2}{A\_{22}} \right] \tag{28}$$

$$E\_h = \frac{1}{t} \left[ A\_{11} - \frac{A\_{12}^2}{A\_{11}} \right] \tag{29}$$

### *2.5. Theoretical Calculation of Design Parameters*

A drive shaft can be idealized as a simply supported beam. To determine the bending natural frequencies for a composite tube, Euler's equation for the lateral vibration of beams can be implemented [57,58]:

$$
\omega\_n = \mathcal{C}\_n \sqrt{\frac{E\_x I}{mL^4}} \tag{30}
$$

or

$$f\_n = \frac{C\_n}{2\pi} \sqrt{\frac{E\_x I}{mL^4}}\tag{31}$$

where *ω<sup>n</sup>* (rad/s) is the angular velocity, *fn* (Hz) is the fundamental (rotational) frequency, *m* is the mass of a drive shaft per unit length, *L* is the length of a drive shaft, and *I* is the second moment of inertia given for a thin-walled tube as:

$$I\_x = \frac{\pi}{4} (r\_0^4 - r\_i^4) \tag{32}$$

where *r*<sup>0</sup> is the outer radius, *ri* is the inner radius, (*ExI*) represents the bending stiffness for a composite tube, and the number *Cn* depends on the boundary conditions of the drive shaft and the number of natural frequency *n*. Numerical values of *Cn* for simply supported boundary conditions [59] are (nπ) 2.

The vibration of parts and assemblies increases as operating conditions become more demanding. Critical speed is defined as the angular rotating shaft speed that is equal to the lowest frequency of its natural vibration. Equations (30) and (31) are used to calculate the natural shaft oscillation (critical speeds) [60].

The main loading of a drive shaft is torsion. Instability may occur at a certain amount of loading. This parameter is critical in the design of drive shafts, and it is important to have an order of magnitude solution before starting the designing [61]. There are various empirical equations in the literature that are based on experimental studies and are used to calculate the torsional buckling load of long thin-walled shafts. The expression of the critical buckling torque for thin-walled orthotropic tubes is given as in [56]:

$$T\_{cr} = (2\pi r^2 t) (0.272) \left[ E\_x E\_h^3 \right]^{1/4} \left( \frac{t}{r} \right)^{3/2} \tag{33}$$

where *r* and *t* are the mean radius and overall thickness of the composite drive shaft, respectively.

### **3. Finite Element Modeling**

In this study, the numerical simulation of a hybrid drive shaft was performed using FEM. The composite drive shaft is considered as a thin-walled tube, balanced, while the stress–strain relationship is linear and elastic. The analysis consists of the following four stages.

### *3.1. Defining Material Properties*

In the current study, epoxy resin, multi-walled CNTs, and carbon/glass fiber were chosen as matrix material, nano material, and reinforcing material, respectively. Material properties of the components constituting the composite material are shown in Table 1.


**Table 1.** Materials properties of the components.

The material properties of CNTs-reinforced unidirectional composite laminas were calculated with two different approaches. In the first approach, concerning e-glass fiber/epoxy resin/CNTs, Equations (7)–(12) were used to calculate its orthotropic mechanical properties. In the second approach, concerning carbon fiber/epoxy resin/CNTs laminas, Equations (7), (11) and (13) were used to calculate the corresponding properties. Additionally, volume fractions of fiber and matrix were considered as 57.65% and 42.35%, respectively. In this way, the composite lamina properties without CNTs inclusions were close to the values referred in [56,58]. Following the approach presented in the previous section, the engineering constants are shown in Table 2. The CNTs volume fractions were considered as 0.00, 0.05, 0.10, 0.20, 0.35, and 0.50.

**Table 2.** Material properties of unidirectional composite laminas.


Elastic constants *E*1, *E*2, *G*12, and *G*<sup>23</sup> were found to increase with the volume fraction of CNTs in the polymer matrix. Specifically, in the first approach, 35%, 582%, 575%, and 607% increases in *E*1, *E*2, *G*12, and *G*23, respectively, were observed for *Vcnt* = 0.50 compared with the laminated composite drive shaft without CNTs inclusions. In the second approach, 9%, 174%, 222%, and 182% corresponding increases were observed. In the first approach, the Poisson's ratio ν<sup>12</sup> seemed to increase as *Vcnt* increased up to 0.35 and then remain stable, while in the second approach, it increased as the *Vcnt* increased up to 0.20 vol% and then remained stable. Finally, in the first approach, the Poisson's ratio ν<sup>23</sup> increased as the *Vcnt* increased up to 0.10 and then decreased, while, in the second approach, it decreased with *V*cnt.

### *3.2. Configuration and Finite Element Type*

A composite drive shaft was considered with a length of 1730 mm and mean radius of 50.8 mm, which consists of four layers, stacked up as a [±45◦ glass/0◦ carbon/90◦ glass] layup. The selected thicknesses for each layer were: +45◦ glass = 0.1905 mm, −45◦ glass = 0.1905 mm, 0◦ carbon = 0.635 mm, 90◦ glass = 1.016 mm.

The finite element modeling was implemented in ABAQUS CAE code. The drive shaft was discretized by using the S4R element type. This is a 3D 4-node, quadrilateral, stress/displacement shell element with reduced integration and a large-strain formulation. These elements allow transverse shear deformation and account for finite membrane strains and arbitrarily large rotations. Reduced integration usually provides more accurate results and significantly reduces the running time, especially in three dimensions. The total number of elements with a global element size *l <sup>e</sup>* = 10 mm was calculated equal to 5536 elements.

### *3.3. Boundary Conditions*

Appropriate boundary conditions must be applied to the drive shaft finite element model according to the analysis type, i.e., modal or elastic buckling analysis. In the modal analysis, boundary conditions are considered to be pinned-pinned, i.e., the drive shaft is simply supported in order to obtain the mode shapes and the corresponding natural frequencies.

In the buckling analysis, the model is fully restrained (fixed) at one end, while the drive shaft is subjected to a torsion load (moment) at the other end. In this way, the buckling modes and loads can be evaluated using a linear buckling analysis. The lowest value is the critical torsional buckling load.

### *3.4. Meshing Sensitivity Analysis*

The effect of global size *l <sup>e</sup>* on the natural frequency and buckling torque is presented in Figure 2a,b, respectively. It is obvious that a relatively coarse mesh was enough to reach the first 10 natural frequencies with reasonable accuracy, while a denser mesh was required in order for the model to converge in a specific magnitude of the critical buckling torque. Specifically, in Figure 2b, the solution seemed to converge with a global element size equal to 10 mm. This is the reason that, here, the global size *l <sup>e</sup>* of the element in the finite element analysis was chosen to be *l <sup>e</sup>* = 10 mm.

**Figure 2.** Mesh sensitivity analysis diagram of a hybrid composite driveshaft consisting of 4 layers, stacked up as a [±45◦ glass/0◦ carbon/90◦ glass] layup, computed for lg = 1730 mm, r = 50.8 mm, and w.thk = 2.032 mm, for (**a**) modal analysis, and (**b**) linear bucking analysis, where the control parameter *l <sup>e</sup>* was increased. To detect the appropriate global size, the model was discretized by using the S4R type of element, with the total No of 5536 elements used for *l <sup>e</sup>* = 10 mm.

### **4. Results**

### *4.1. Validation*

To evaluate the validity of the proposed method in terms of natural frequency and torsional buckling load, a comparison with the analytical and numerical results available in [56] was performed concerning the drive shaft of the configuration presented in Section 3.2. The comparison is presented in Table 3 and shows a reasonable agreement between the results of the different methods. Note that the present model used engineering constants computed by the approach presented in the previous section that slightly differ with the corresponding values of [58]. Nevertheless, the proposed finite element model has also been tested using identical properties used in [58], and the deviations on the results were calculated to be lower than 1%.



Figure 3 illustrates a set of eight natural frequencies and their corresponding mode shapes. The results also agree with the results presented in [58] in terms of both mode shapes and natural frequencies values. The greatest discrepancies between the two methods were observed in the breathing vibration modes; however, these modes were not involved in the basic (first three) mode shapes of vibration of the specific configuration of the composite shaft.

Furthermore, the fundamental natural frequency of a drive shaft consisting of four layers stacked as [±45◦ glass/0◦ carbon/90◦ glass] for various fiber orientations was evaluated. The variation is presented in Figure 4. The results were also compared with the results of the study [58] presenting, again, a reasonable agreement. It is observed that the fundamental frequency (Figure 4a) increased as the fiber angle of the first two ±θ◦ glass layers tended to 0◦. However, the optimum configuration may be determined taking into consideration the other design parameters as well. In Figure 4b, it is observed that the fundamental frequency decreased, as the fiber angle of the third 0◦ carbon layer tended to 90◦. The results indicate that the drive shaft lost 55% of its natural frequency when the carbon fibers oriented at 90◦ in the hoop direction instead of 0◦. In order to increase the modulus of elasticity in the longitudinal direction of the drive shaft, the carbon fibers layer, with its high modulus, must be oriented at zero angle.

**Figure 3.** Set of the first eight mode shapes of natural frequency of a simply supported composite drive shaft. To detect deformation and natural frequencies, FEM analysis was conducted in Abaqus CAE. The corresponding frequency values were as follows: (**a**) 91.65, (**b**) 326.88, (**c**) 635.73, (**d**) 442.34, (**e**) 334.69, (**f**) 334.91, (**g**) 402.28, and (**h**) 1485.60 Hz.

**Figure 4.** Fundamental frequency as a function of fiber orientation angle (θ◦). A hybrid 372 composite drive shaft having four layers stacking of (**a**) [±θ◦glass/0◦carbon/90◦glass] and (**b**) [±45◦glass/θ◦carbon/90◦glass] was analyzed for length = 1730 mm, radius = 50.8 mm, and wall thickness = 2.032 mm, for modal analysis, where the fiber orientation and the control parameter θ◦ [deg] were transformed. Figures (**a**,**b**) show the effect of fiber orientation angle on f [Hz] by changing the first two layers and the third (carbon fiber) layer, respectively. To ensure the validity of the process, identical material properties were employed. The present approach and FEM [56] are presented by blue solid lines and dashed red lines, respectively.

### *4.2. Effect of CNTs Inclusion on Natural Frequencies and Critical Speed of the Composite Drive Shaft*

In order to determine the effect of CNTs on the composite drive shaft on natural frequency, a parametric finite element analysis was conducted. The results are depicted in Figure 5.

**Figure 5.** Natural frequencies as a function CNTs volume fraction. A hybrid composite drive shaft of stacking [±45◦ glass/0◦ carbon/90◦ glass] was analyzed for length = 1730 mm, radius = 50.8 mm, and wall thickness = 2.032 mm, for modal analysis, where the volumetric concentration of nanoparticles was increased. The gradient of the lines shows the effect of CNTs volume fraction on the natural frequencies computed for the following ratios: 0, 0.05, 0.10, 0.20, 0.35, and 0.50.

It is observed that as the volume fraction of CNTs increased, the natural frequency tended to increase. Specifically, the enhancement of the values of this set of natural frequencies reached 59%, 69%, 80%, 133%, 92%, 94%, 89%, and 55% for the 1st bending,

2nd bending, 3rd bending, 1st torsional, 1st breathing, 2nd breathing, 3rd breathing, and 1st axial natural frequency, respectively, concerning a 0.50 volume fraction of CNTs.

Figure 6 shows the variation in the critical speed, for different fiber orientation angles, and for various CNT volume fractions in the epoxy resin matrix. For volume fractions of CNTs equal to 0.50 into the hybrid composite drive shaft, the critical speed presented a significant increase in its value, approximately 60%, as depicted in Figure 6a,b, respectively.

**Figure 6.** Critical speed as a function of fiber orientation angle (θ◦) for various CNT volume fractions. A hybrid composite drive shaft having four layers stacking of (**a**) [±θ◦ glass/0◦ carbon/90◦ glass] and (**b**) [±45◦ glass/θ◦ carbon/90◦ glass] was analyzed for length = 1730 mm, radius = 50.8 mm, and wall thickness = 2.032 mm, where the fiber orientation and volumetric concentration of nanoparticles were changed. Figures (**a**,**b**) show the effect of fiber orientation angle on the critical speed by changing the first two layers and the third (carbon fiber) layer, respectively. The curves show the effect of CNTs volume fractions on the critical speed, computed for CNT volume fractions: 0, 0.05, 0.10, 0.20, 0.35, and 0.50.

### *4.3. Effect of CNTs Inclusion on Buckling Torque*

The first mode shape of buckling torque, as derived from the finite element analysis, is presented in Figure 7. The first buckling mode was similar to the corresponding one predicted in [61]. In this analysis, the computed load was close to the collapse load. The output from the analysis is a factor (eigenvalue) that is multiplied by the actual magnitude of the applied load in order to compute an estimation of the critical torque. To determine the variation in the composite drive shaft under buckling torque due to the presence of CNTs into the composite shaft, a parametric finite element investigation was also performed. The results are presented in Figure 8. It was observed that the performance of the critical buckling torque significantly improved with the amount of added CNTs in the matrix of the composite drive shaft.

This may be expected due to the increase in the mechanical properties of the composite lamina as the CNTs volume fraction is increased. Moreover, the variation in the critical buckling torque, for different fiber orientation angles, is also depicted. With the addition of CNT into the pure composite drive shaft, the critical buckling torque showed an enhancement in its value. Specifically, in Figure 8a, the critical buckling torque of the lamination [±0, 0, 90] increased up to 27%, 50%, 89%, 140%, and 187% for *Vcnt* = 0.05, 0.10, 0.20, 0.35, and 0.50, respectively. The critical buckling torque of the lamination [±45, 0, 90] increased up to 19%, 35%, 66%, 107%, and 145% for *Vcnt* = 0.05, 0.10, 0.20, 0.35, and 0.50, respectively. The critical buckling torque of the lamination [±90, 0, 90] increased 12%, 21%, 37%, 58%, and 79% for *Vcnt* = 0.05, 0.10, 0.20, 0.35, and 0.50, respectively. Moreover, it seems that the critical buckling torque was not considerably affected by the fiber orientation angle, as the amount of added CNTs increased. Finally, in Figure 8b, the critical buckling torque

of the lamination [±45, 0, 90] increased up to 19%, 35%, 66%, 107% and 145% for *Vcnt* = 0.05, 0.10, 0.20, 0.35, and 0.50, respectively. Additionally, we notice that the critical buckling torque increased as the fiber angle of the third ply (carbon fiber) tended to 0◦ or 90◦. On the contrary, the lowest value was observed at 45◦.

**Figure 7.** First buckling mode shape of a hybrid composite driveshaft, stacked as [±45◦ glass/0◦ carbon/90◦ glass], computed for length = 1730 mm, radius = 50.8 mm, and wall thickness = 2.032 mm. To detect the deformation and critical buckling load, eigenvalue buckling analysis was conducted, where the model was fully fixed at one end and subjected to torsion load at the other end. The corresponding buckling torque value from Equation (33) was 2149 N·m and from FEA was 1730 N·m.

**Figure 8.** Critical buckling torque as a function of fiber orientation angle (θ◦) for various 433 CNT volume fractions. A hybrid composite drive shaft having four layers stacking of (**a**) [±θ◦ glass/0◦ carbon/90◦ glass] and (**b**) [±45◦ glass/θ◦ carbon/90◦ glass] was analyzed with length = 1730 mm, radius = 50.8 mm, and wall thickness = 2.032 mm. To obtain the critical buckling torque eigenvalue, buckling analysis was conducted.

### **5. Conclusions**

In this study, the effect of MWCNTs on a hybrid fiber-reinforced composite automotive drive shaft was investigated. The modified epoxy resin matrix with different CNTs volume fractions was theoretically modeled. The CNTs-modified matrix presented enhanced material properties compared to the pure epoxy matrix system, affecting the overall behavior of the drive shaft. According to the results, the following findings can be noticed:

• The elastic constants *E*1, *E*2, *G*12, and *G*<sup>23</sup> could be increased up to 35%, 582%, 575%, and 607% (e-glass-reinforced polymer lamina) and 9%, 174%, 222%, and 182% (carbon

fiber-reinforced polymer lamina), respectively, due to the 0.50 volume fraction of CNTs in the polymer matrix for the specific drive shaft configuration considering uniform dispersion.


Eventually, the enhanced characteristics of CNTs-based nanocomposites have the potential to be employed for the design of composite structures in the areas of automotive, aerospace, construction, military, etc., in terms of improving stiffness or strength, in addition to its low weight. This computational study did not investigate the effects of CNT dispersion and distribution into the polymer matrix. Issues such as agglomeration are going to be modeled/studied theoretically and/or experimentally in future work. The obtained results can be considered as an upper bound limit of the potential enhanced mechanical performance. It seems that CNTs can play an important role as a significant secondary filler and the proposed approach may be used as a design tool.

**Author Contributions:** Conceptualization, S.K.G.; methodology, S.K.G.; software, P.A.A.; validation, P.A.A. and S.I.M.; formal analysis, S.K.G.; investigation, P.A.A.; resources, P.A.A.; data curation, P.A.A.; writing—original draft preparation, S.K.G.; writing—review and editing, S.K.G.; visualization, S.I.M.; supervision, S.K.G.; project administration, S.K.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

