A Modified Couple Stress Elasticity for Non-Uniform Composite Laminated Beams Based on the Ritz Formulation

Due to the large application of tapered beams in smart devices, such as scanning tunneling microscopes (STM), nano/micro electromechanical systems (NEMS/MEMS), atomic force microscopes (AFM), as well as in military aircraft applications, this study deals with the vibration behavior of laminated composite non-uniform nanobeams subjected to different boundary conditions. The micro-structural size-dependent free vibration response of composite laminated Euler–Bernoulli beams is here analyzed based on a modified couple stress elasticity, which accounts for the presence of a length scale parameter. The governing equations and boundary conditions of the problem are developed using the Hamilton’s principle, and solved by means of the Rayleigh–Ritz method. The accuracy and stability of the proposed formulation is checked through a convergence and comparative study with respect to the available literature. A large parametric study is conducted to investigate the effect of the length-scale parameter, non-uniformity parameter, size dimension and boundary conditions on the natural frequencies of laminated composite tapered beams, as useful for design and optimization purposes of small-scale devices, due to their structural tailoring capabilities, damage tolerance, and their potential for creating reduction in weight.


Introduction
In the last decades, composite structures and materials have received an increased interest in many industries such as the aerospace, automotive, biomedical, architectural, mechanical, and civil sectors [1], due to their high mechanical performances. In particular, micro/nano-scale mechanical structures usually feature a characteristic size of micron or submicron order, e.g., micro/nano-beams, and micro-nano-cylinders, largely used in micro-and nano-electromechanical devices (MEMS and NEMS). Several experimental evidences in literature, have revealed that the behavior of micro-structures is size-dependent [2][3][4][5]. Thus, a large number of works has been recently published to conceive novel structural solutions, systems, and devices, while adopting different types of reinforcement phase, such as graphene nanoplatelets [6][7][8][9][10][11][12][13][14], or carbon nanotubes [15][16][17][18][19]. Among a large variety of numerical strategies, higher order theories represent the most useful tool for the investigation of the static and dynamic response of materials at different scales [20][21][22][23][24]. Classical theories, indeed, similar CST-based assumptions as in Ref. [35], which are here generalized to handle composite laminated Euler-Bernoulli beams with a more complicated tapered geometry and different boundary conditions. Among different numerical approaches, in the present work we apply a Ritz-type solution with harmonic trial functions to solve the problem, whose stability and accuracy is verified through a systematic investigation. In line with predictions from the literature [65][66][67][68][69][70][71][72][73], the Rayleigh-Ritz method, represents an efficient tool for the analysis of the structural behavior of beams, whose accuracy and stability are well known to be related to the selected trial functions. The trial functions must satisfy the enforced boundary conditions. When this condition is not fulfilled, the Lagrangian multipliers and penalty method could be adopted to handle arbitrary boundary conditions. This approach, however, can cause an overall increase in dimension for both the stiffness and mass matrices, with a consecutive increase in the computational cost. Therefore, in the present work we first check for the stability of the numerical solution for the selected harmonic trial functions, by means of a systematic investigation. The numerical study also aims at evaluating the sensitivity of the response to different geometrical and/or mechanical parameters, which could be of great interest for design purposes in practical engineering application, and could serve for future studies on nonuniform beams and devices.
The outline of the paper is as follows: in Section 2 we introduce the mathematical problem for tapered nanobeams, which is solved numerically by means of the Ritz method in Section 3. The numerical examples and applications are discussed comparatively in Section 4 for different mechanical and geometrical parameters. Finally, in Section 5, we draw the main conclusions of our work.

Theory and Mathematical Problem
Let us consider the orthotropic non-uniform nanobeam in Figure 1, with length , constant thickness ℎ, variable width ( ), in a Cartesian coordinate system (x,y,z). Based on the Timoshenko beam theory, the displacement field u , is defined by its components, as follows [35,61]: Based on the Timoshenko beam theory, the displacement field u, is defined by its components, as follows [35,61]: where u 0 (x, t) and w 0 (x, t) are the axial and transverse displacements of an arbitrary point of the mid-plane along the x-and z-directions, respectively, whereas φ(x, t) is the angle of rotation around the y-axis of the cross section, that will be defined as φ(x, t) = ∂w 0 (x, t)/∂x for an Euler-Bernoulli formulation. Based on the NMCST, the rotational field θ = 1/2curlu is defined by the following components: The non-null components of the strain tensor ε for the kth ply of a laminated beam, are governed by the following kinematic relations: where γ k xz becomes equal to zero according to Euler-Bernoulli theory. In line with the NMCST proposed in [35], we introduce the constitutive relations for the kth ply of a laminated micro-composite beam, in the global system of coordinates, where two length scale parameters are introduced, l 2 kb and l 2 km for fibers and matrix in the kth lamina, respectively. More specifically, l 2 kb refers to the micro-scale material constant of an arbitrary fiber rotating in the y − z plane, where the fiber is considered as the impurity affecting the rotational equilibrium; l 2 km stands for the micro-scale material constant within the matrix rotating about the impurity in the x − z plane.
Thus, the stress-strain relations in the global coordinate system, are expressed in compact form as: where σ k = σ k x σ k y τ k xz τ k yz m k xy m k yx T , ε = ε x ε y γ xz γ yz χ xy χ yx T , m ij stand for the modified couple stresses, and Q k = T kT C k T k depends on the coordinate transformation matrix T k and on the elastic properties matrix C k , defined as follows: In the matrix (5), m = cos ϕ k , n = sin ϕ k ,ϕ k is the fiber angle of a layer with respect to the x-axis, while the elastic stiffness components C ij in matrix (6) are defined as in Reference [35]: with G ij and E i the shear and normal elastic modulus, respectively, and ν ij the Poisson ratios. Once the coordinate transformation from a local to the global system is performed, the constitutive relations take the following form: where the elastic coefficients are defined as [35]: Starting with the above-mentioned constitutive relations for composite laminated beams based on the NMCST, we determine the governing equations of motion by means of the Hamilton's principle. In absence of external forces acting on the structure, the total potential energy Π takes the following form: U and K being the strain energy and the kinetic energy, respectively. More specifically, the strain energy of the beam is defined, in the domain V, as follows: which is combined to Equations (3) and (4) to yield the following expression: Molecules 2020, 25, 1404 Under the Euler-Bernoulli assumption of the type φ(x, t) = ∂w 0 (x, t)/∂x, Equation (12) becomes: where and b(x) refers to the non-uniform width, whose variation is defined as: where b 0 is the width of the tapered beam, at x = 0, and N is the exponential non-uniform parameter. The kinetic energy in Equation (10) is expressed as follows: For a Euler-Bernoulli beam formulation, Equation (16) becomes as follows: with Molecules 2020, 25, 1404 By a combination of Equations (10), (13), (17) we get the following expression for the total energy for the Euler-Bernoulli beam:

The Rayleigh-Ritz Procedure
The Rayleigh-Ritz method, with two different exponential trial functions, is here applied to approximate the displacement field as proposed by Nguyen et al. [72], and determine the solution of the problem. Thus, the kinematic quantities are approximated as follows: where ω is the natural frequency, i 2 = −1 refers to the imaginary unit, u j , w j are the unknowns of the problem, and ψ j are the trial functions which depend on the selected boundary conditions. In the present study we consider two different types of boundary conditions, namely simply supports (S-S) and clamped-free (C-F) supports, such that the following trial functions are assumed [71]: Upon substitution of Equations (20), (21) into Equation (19), and by using the Lagrange's equations, we get the following relation: p j being the values of u j , w j , that describe the vibration response of the tapered beam structure. After some mathematical manipulation, the generalized eigenvalue problem gets the following form where K and M stand for the stiffness and mass matrix, respectively, whose components are defined as follows: The natural frequencies of the orthotropic nanostructure are, finally, determined through the enforcement of the following condition:

Numerical Results and Discussion
In this section, we present the results of different numerical examples, selected to test the accuracy of the formulation with respect to the available literature, and the sensitivity of the free vibration response to the boundary conditions, length-scale parameter, non-uniformity parameter, or size dimension.
For validation purposes, we compute the first five natural frequencies for a S-S three-layer [90 • ,0 • ,90 • ] microbeam, with the following geometrical properties: b = h = 25 × 10 −6 m, L = 200 × 10 −6 m. The mechanical properties of the material are assumed as in Reference [35], i.e., E 2 = 6.9 × 10 9 Pa, E 1 = 25E 2 , G 12 = G 13 = 0.5E 2 , G 23 = 0.2E 2 , ν 12 = ν 13 = ν 23 = 0.25, ρ = 1.578 kg/m 3 . Table 1 summarizes the results for the first five frequency parameters, and different values of m. As clearly visible in Table 1, m = 5 represents a convergence point for the numerical computation of the natural frequencies. This value of m is assumed hereafter for the parametric study. A further comparative example is chosen to assess the capabilities of the present formulation, namely, an isotropic S-S uniform nanobeam, as proposed originally by Chen and Li [35]. The first five natural frequencies computed with our formulation are compared to predictions by Chen and Li [35], as summarized in Table 2, for a different length-scale parameter. A good agreement with the available literature is observed, which confirms the accuracy of the proposed formulation, along with a general increase of each natural frequency for an increased length-scale parameter. Small differences between our predictions and the ones in Reference [35] are noticed for an increased length scale parameter. This is mainly related to the different basic assumptions considered in the two works, namely a Euler-Bernoulli beam model instead of a Timoshenko-based formulation. In agreement with findings by Reference [35], it seems that an Euler-Bernoulli-based formulation gets higher natural frequencies than a Timoshenko-based theory. After the preliminary validation, we perform a parametric analysis of the vibration response for an orthotropic non-uniform nanobeam under two different sets of boundary condition, while including the effects of size, length-scale and non-uniformity. A three layer [90 • ,0 • ,90 • ] non-uniform nanobeam is considered, with h = 10 nm, b 0 = 2 h, different values of L, and material properties stemming from Reference [23], i.e., E 2 = 13.67 GPa, E 1 = 37.41 GPa, G 12 = 6.03 GPa, G 13 = 6.03 GPa, G 23 = 6.67 GPa, ν 12 = ν 13 = ν 23 = 0.3, ρ = 1938.9 kg/m 3 . Tables 3 and 4 illustrate the main results in terms of natural frequency for different size ratios, L/h, non-uniformity parameter, Nh, length scale, l, for a S-S and C-F non-uniform nanobeam, respectively. Based on these tables, an increased length scale and a decreased size ratio leads to an overall increase of the natural frequencies. A non-monotonic behavior, instead, is exhibited by the natural frequencies for an increasing non-uniformity Nh, while keeping fixed the other parameters. In Figure 2, we plot the variation of the first and the fifth natural frequencies for a three-layer [90 • ,0 • ,90 • ] Euler-Bernoulli S-S non-uniform beam with size ratio L/h, and fixed values of h = 10 nm, b 0 = 2 h, Nh = 0.5. The NMCST under the assumption of l = 1 nm is here compared to the classical approach (i.e., for l = 0 nm). As clearly shown in Figure 2, the classical theory predicts lower values of natural frequencies with respect to the NMCST here proposed, whereby, both natural frequencies (ω 1 , ω 5 ) decrease for increased geometrical lengths of nanobeams. The parametric study is thus repeated for a C-F nanobeam, whose results are plotted in Figure 3 in terms of the natural frequencies ω 1 , ω 5 , while assuming the same geometrical and mechanical parameters as in the previous investigation. Additionally, in this case, the NMCST yields higher values of natural frequencies compared to a classical approach. The main difference between the two approaches, in this case, is less pronounced because of the lower deformability of the C-F nanobeam compared to a S-S boundary condition.  Nh=3 or 2, depending on the selected boundary condition, with a drastic reduction up to a null asymptotic value. A monotone increase of the natural frequency is also observed for an increasing length scale parameter l, at least for lower values of ℎ, whose variation is finally visualized in Figure  5a,b, for a S-S-and C-F nanostructure, respectively. Based on the last results, it seems that uniform beams are more sensitive to the length parameter, compared to tapered geometries, which could be accounted for design purposes of nanodevises.   asymptotic value. A monotone increase of the natural frequency is also observed for an increasing length scale parameter l, at least for lower values of ℎ, whose variation is finally visualized in Figure  5a,b, for a S-S-and C-F nanostructure, respectively. Based on the last results, it seems that uniform beams are more sensitive to the length parameter, compared to tapered geometries, which could be accounted for design purposes of nanodevises.   In Figure 4a,b, we plot the first natural frequency vs. the non-uniformity parameter Nh, for the S-S and C-F non-uniform nanobeam, along with different length scale parameters, namely, l = 0; 0.1; 0.5; 1, and a fixed geometry h = 10 nm, b 0 = 2h, L/h = 10. As shown in Figure 4, for both sets of boundary conditions, the natural frequency decreases with an increased non-uniformity parameter (see the zoom-ups of Figure 4). Moreover, the natural frequency decreases monotonically between Nh = 0 and Nh = 3 or 2, depending on the selected boundary condition, with a drastic reduction up to a null asymptotic value. A monotone increase of the natural frequency is also observed for an increasing length scale parameter l, at least for lower values of Nh, whose variation is finally visualized in Figure 5a,b, for a S-S-and C-F nanostructure, respectively. Based on the last results, it seems that uniform beams are more sensitive to the length parameter, compared to tapered geometries, which could be accounted for design purposes of nanodevises.

Conclusions
In this work, we employ a novel modified couple stress theory for studying the vibration response of laminated composite non-uniform beams under two different boundary conditions. The problem is tackled with the Rayleigh-Ritz formulation, here proposed as promising numerical approach to predict the size-dependent responses of micro composite beams. This is verified through a comparative study with the available literature, at least for uniform geometries. A parametric investigation is, thus, repeated systematically to check for the sensitivity of the vibration behavior for non-uniform nanobeams with different geometrical shape, non-uniformity parameter, length scale parameter, and boundary condition. The numerical outcomes show that the length-to-thickness ratio, non-uniformity, and boundary condition, play a key role in the vibration response of the nanostructure, compared to the length-scale sensitivity. More specifically, an increased length scale and a decreased size ratio yields an overall increase in the natural frequencies, along with an increased stiffness. As expected, a classical theory predicts lower values of natural frequencies with respect to the NMCST here proposed, whereby, an increased geometrical length of the nanobeams yields an overall decrease in the natural frequencies and structural stiffness. In addition, an increased non-uniformity in the beam gets lower natural frequencies. This means that the non-uniformity parameter of a tapered beam could enable a tailorable stiffness and vibration response, depending on the design requirements. These conclusions could be of interest for the nanotechnology community, as well as for design purposes and optimization processes of many engineering nanodevices, nanoelectronics, or nanosensors. The basic notions of the formulation here proposed, could be also used to treat other mechanical aspects, such as buckling problems or fracture mechanics problems of tapered beams.