Buckling Behavior of Nanobeams Placed in Electromagnetic Field Using Shifted Chebyshev Polynomials-Based Rayleigh-Ritz Method

In the present investigation, the buckling behavior of Euler–Bernoulli nanobeam, which is placed in an electro-magnetic field, is investigated in the framework of Eringen’s nonlocal theory. Critical buckling load for all the classical boundary conditions such as “Pined–Pined (P-P), Clamped–Pined (C-P), Clamped–Clamped (C-C), and Clamped-Free (C-F)” are obtained using shifted Chebyshev polynomials-based Rayleigh-Ritz method. The main advantage of the shifted Chebyshev polynomials is that it does not make the system ill-conditioning with the higher number of terms in the approximation due to the orthogonality of the functions. Validation and convergence studies of the model have been carried out for different cases. Also, a closed-form solution has been obtained for the “Pined–Pined (P-P)” boundary condition using Navier’s technique, and the numerical results obtained for the “Pined–Pined (P-P)” boundary condition are validated with a closed-form solution. Further, the effects of various scaling parameters on the critical buckling load have been explored, and new results are presented as Figures and Tables. Finally, buckling mode shapes are also plotted to show the sensitiveness of the critical buckling load.


Introduction
In recent decades, with the advent of technological advancement, small-scale structures like nanoclock, nano-oscillator, nanosensor, NEMS (Nano-Electro-Mechanica System), and so forth have found tremendous attention due to their various applications. In this scenario, the study of dynamical behaviors of nanostructures is important and a need of the hour. Due to the small size of nanostructures, experimental analysis of such structures is always challenging and difficult as it requires enormous experimental efforts. Moreover, classical mechanics or continuum mechanics fail to address the nonlocal effect also. In this regard, nonlocal continuum theories were recently found to be useful for the modeling of micro-and nanosized structures. Various researchers developed different nonlocal or nonclassical continuum theories to assimilate the small-scale effect, such as strain gradient theory [1], couple stress theory [2], micropolar theory [3], nonlocal elasticity theory [4], and so on. Out of all these theories, nonlocal elasticity theory of Eringen has been broadly used for the dynamic analysis of nanostructures. Few studies regarding the vibration and buckling of beam, membrane, and nanostructures such as nanobeam, nanotube, nanoribbon, and so forth can be found in [5][6][7][8][9][10][11][12][13][14].
Wang et al. [15] studied buckling behavior of micro-and nanorods/tubes with the help of Timoshenko beam theory, where small-scale effect was addressed by using the nonlocal elasticity theory of Eringen. Emam [16] incorporated nonlocal elasticity theory to analyze the buckling and

Proposed Model for Electromagnetic Nanobeam
In this study, the nanobeam with length L and diameter d is placed in an electromagnetic field with the electric field intensity as E and the magnetic flux density as B. The schematic diagram for continuum model of the nanobeam is shown in Figure 1. Then, by Ohm's law, the current density (J) of the system due to the induced current (because of Lorentz force) is given as [31] where σ 0 is the electrical conductivity, µ 0 is the magnetic permeability of free space, and H is the magnetic field strength. By neglecting the electric field intensity, the nanobeam experiences a magnetic force or Pondermotive force which is denoted by f em and can be expressed as [31] f em = J × B = σ 0 (E +w 0 × µ 0 H) × µ 0 H = σ 0 µ 2 0 H 2 w 0 (2)  According to Euler-Bernoulli beam theory, the displacement field at any point may be stated as [32] ( Here, 1 u and 3 u represent displacements along x and z directions, respectively, and where xx σ is the normal stress, xx ε is the normal strain, and  According to Euler-Bernoulli beam theory, the displacement field at any point may be stated as [32] Here, u 1 and u 3 represent displacements along x and z directions, respectively, and w 0 (x, t) denotes the transverse displacement of the point on the mid-plane of the beam. The strain-displacement relation may be expressed as Under the framework of Euler-Bernoulli nanobeam, the variation of strain energy (δU) and the variation of work done by external force (δW e ) are presented as where σ xx is the normal stress, ε xx is the normal strain, and M xx = A σ xx z dA is the bending moment of nanobeam. The Hamilton's principle for the conservative system is stated as Nanomaterials 2019, 9, 1326 4 of 15 Substituting Equations (5)-(7) and setting δ = 0, we have The equation of motion for buckling behavior can be obtained as For an isotropic nonlocal beam, the nonlocal elasticity theory of Eringen can be expressed as [4] 1 − µ ∂ 2 ∂x 2 σ xx = Eε xx (10) where µ = (e 0 a) 2 is the nonlocal parameter, E is Young's modulus. Here e 0 and a denote material constant and internal characteristic length, respectively. Multiplying Equation (10) by zdA and integrating over A, the nonlocal constitutive relation for Euler-Bernoulli nanobeam may be expressed as where I = A z 2 dA, is the second moment of area. Using Equation (9) in Equation (11) and rearranging, the nonlocal bending moment can be obtained as Equating the nonlocal strain energy with work done by an external force, we obtain Substituting Equation (12) in Equation (9), we obtain the governing equation of motion as Let us define the following nondimensional parameters Incorporating the above nondimensional parameters in Equations (13) and (14), we have

Shifted Chebyshev Polynomials-Based Rayleigh-Ritz Method
Chebyshev polynomials of the first kind (C n (X)) are a sequence of orthogonal polynomials with X ∈ [−1 1], and few terms of the sequence are defined as In order to solve Equation (12), Rayleigh-Ritz method is implemented along with Chebyshev polynomials of the first kind as shape function. For more details about the Rayleigh-Ritz method, one may refer to the books [33,34]. The main advantages of using Chebyshev polynomials over algebraic polynomials 1, X, X 2 , X 3 , . . . X n are the orthogonal properties of Chebyshev polynomials, which reduce the computational effort, and for the higher value of n (n > 10), the system avoids ill-conditioning. Since the domain of the nanobeam lies in [0 1], the Chebyshev polynomials must be reduced to shifted Chebyshev polynomials of the first kind (C * n (X)) with X ∈ [0 1]. This is achieved by transforming X → 2X − 1 , and there exists a one-to-one correspondence between [0 1] and [−1 1]. Accordingly, the first few terms of shifted Chebyshev polynomials of the first kind (C * n (X)) can be written as, The transverse displacement function (W(X)) as per the Rayleigh-Ritz method can now be expressed as where a i s are unknowns, C * i−1 are the shifted Chebyshev polynomials of the index i − 1, N is the number of terms required to obtain the result with the anticipated accuracy, p and q are the exponents which decide the boundary conditions, as given in Table 1. Table 1. Values of p and q for different boundary conditions [33,34].

Boundary Conditions
Replacing Equation (19) into Equation (15), and minimizing the buckling load parameter with respect to the coefficients of the admissible functions (i.e., a i s, i = 1, 2, 3 . . . N ), we obtain the generalized eigenvalue problem as

Closed-Form Solution for P-P Boundary Condition Using Navier's Technique
Navier's technique has been employed to find a closed-form solution for the Pined-Pined (P-P) boundary condition. As per the Navier's technique, the transverse displacement (W) may be expressed as [23][24][25] In which W n , and ω n are the displacement and frequency of the beam. Now, by substituting Equation (21) in Equation (16), the buckling loadP for Pined-Pined (P-P) boundary condition is calculated asP

Numerical Results and Discussion
Shifted Chebyshev polynomials-based Rayleigh-Ritz method has been employed to convert Equation (15) into the generalized eigenvalue problem given in Equation (20). MATLAB codes have been utilized to solve the generalized eigenvalue problem and to compute the critical buckling load parameter. Likewise, Navier's technique has been adopted to find a closed-form solution for the P-P boundary condition, which is demonstrated in Equation (22). In this regard, the following parameters are taken from [15] for computation purpose E = 1 Tpa, d = 1 nm, and L = 10 nm.

Validation
Results of the present model were authenticated by two ways, firstly, by matching with the numerical results given by Wang et al. (2006) for "P-P, C-P, C-C, and C-F" boundary conditions and secondly, Pined-Pined results were compared with the closed-form solution obtained by Navier's technique. For this purpose, the Hartmann parameter (H a ) in the present model was taken as zero, and the critical buckling load parameters (P cr ) for "P-P, C-P, C-C, and C-F" boundary conditions were taken into investigation. Comparisons of critical buckling load (P cr ) are presented in Table 2. Similarly, Table 3 illustrates the comparison of the P-P boundary condition with the Navier's results, with E = 1 TPa, d = 1 nm, L = 10 nm, and H a = 1. From these comparisons, it is evident that the critical buckling loads of the present model are on a par with [15] in the particular case and with Navier's solution for the P-P boundary condition. Table 2. Comparison of "Critical buckling load" (P cr ) in nN with [15] for L d = 10.

(a) Comparison of P-P and C-P boundary conditions
e 0 a P-P C-P

Convergence
A convergence study has been performed to know the number of terms needed to obtain the results of critical buckling load parameters (P cr ) and verify the present model using the Rayleigh-Ritz method. In this regard, H a = 1, L = 10, and e 0 a = 1 were taken for computation purpose. Both the tabular and graphical results were noted for "P-P, C-P, C-C, and C-F" boundary conditions, which are demonstrated in Table 4 and Figure 2, respectively. The C-F boundary condition is converging faster with N = 5, whereas other edges such as P-P, C-P, and C-C take N = 7 for acquiring the desired accuracy. These results revealed that both the model and the results are useful regarding the present investigation. Table 4. Effect of no. of terms (N) on critical buckling load (P cr ) with H a = 1, L = 10, and e 0 a = 1.
N P-P C-P C-C C-F

Influence of Small Scale Parameter
This subsection is dedicated to investigating the influence of a small scale parameter ( ) a e 0 on critical buckling load parameters and the critical buckling load ratio. The four frequently used boundary conditions such as "P-P, C-P, C-C, and C-F" were taken into consideration with 7 = N , 10 = L , and 2 = a H . Tabular and graphical results are illustrated in Table 4 and Figures 2 and 3 for different a e 0 (0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5). Table 5  qualitatively on buckling load. From Table 5 and Figure 3, it is observed that critical buckling load is decreasing with an increase in small scale parameter ( ) a e 0 , and this decline is more

Influence of Small Scale Parameter
This subsection is dedicated to investigating the influence of a small scale parameter (e 0 a) on critical buckling load parameters and the critical buckling load ratio. The four frequently used boundary conditions such as "P-P, C-P, C-C, and C-F" were taken into consideration with N = 7, L = 10, and H a = 2. Tabular and graphical results are illustrated in Table 4 and Figures 2 and 3 for different e 0 a (0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5). Table 5 and Figure 3 represent the variation of small scale parameter (e 0 a) on critical buckling load, whereas Figure 4 demonstrates the variation of small scale parameter (e 0 a) on the critical buckling load ratio. The critical buckling load ratio may be defined as the ratio of critical buckling load calculated using nonlocal theory and classical theory (e 0 a = 0). This critical buckling load ratio acts as an index to estimate the influence of the small scale parameter (e 0 a) qualitatively on buckling load. From Table 5 and Figure 3, it is observed that critical buckling load is decreasing with an increase in small scale parameter (e 0 a), and this decline is more in case of the C-C boundary condition. From Figure 4, it may also be noted that the influence of the small scale parameter is comparatively more in C-C edge and less in C-F edge. Table 5. Effect of small scale parameter (e 0 a) on critical buckling load (P cr ) in nN with N = 7, L = 10, and H a = 2. e 0 a P-P C-P C-C C-F

Influence of Aspect Ratio
The objective of this subsection is to study the impact of aspect ratio ( ) d L on the critical buckling load ( ) cr P with "P-P, C-P, C-C, and C-F" boundary conditions for different d L (5,10,15,20,25,30,35,40,45,50). The effect of aspect ratio has been reported in Table 6 and Figure 5 for , and 2 = a H , which are respectively. From this study, it is essential to note that the critical buckling load decreases with an increase in aspect ratio ( ) d L . This decrease is more consequential for the lower value of aspect ratio.

Influence of Aspect Ratio
The objective of this subsection is to study the impact of aspect ratio (L/d) on the critical buckling load (P cr ) with "P-P, C-P, C-C, and C-F" boundary conditions for different L/d (5,10,15,20,25,30,35,40,45,50). The effect of aspect ratio has been reported in Table 6 and Figure 5 for N = 7, e 0 a = 1, and H a = 2, which are respectively. From this study, it is essential to note that the critical buckling load decreases with an increase in aspect ratio (L/d). This decrease is more consequential for the lower value of aspect ratio. Table 6. Effect of aspect ratio (L/d) on critical buckling load (P cr ) in nN with N = 7, e 0 a = 1, and H a = 2.
L/d P-P C-P C-C C-F    Table 7 and Figure 6 represent the results for the variation of critical buckling load with

Influence of Hartmann Parameter
Hartmann parameter ( ) a H for "P-P, C-P, C-C, and C-F" boundary conditions. From these results, we may note that the critical buckling load decreases with increase in Hartmann parameter, but this drop in critical buckling load is very slow.

Influence of Hartmann Parameter
For the designing of electromagnetic devices, proper knowledge about the effect of electric and magnetic fields on critical buckling load is necessary as it greatly influences the lifespan of electromagnetic devices. In this regard, the effect of Hartmann parameter (H a ) on the critical buckling load (P cr ) has been studied in this subsection for different values of H a (0, 1, 1.5, 2, 2.5, 3, 3.5). Table 7 and Figure 6 represent the results for the variation of critical buckling load with Hartmann parameter (H a ) for "P-P, C-P, C-C, and C-F" boundary conditions. From these results, we may note that the critical buckling load decreases with increase in Hartmann parameter, but this drop in critical buckling load is very slow.

Buckling Mode Shape
Buckling is the state of instability of structures that leads to structural failure. In this circumstance, the buckling mode shape has a vital role in predicting the instability. In this regard, buckling mode shapes were plotted with for different a e 0 (0.5, 1, 1.5, 2). Figures 7-10 show the buckling mode shapes for "P-P, C-P, C-C, and C-F" boundary conditions, respectively. From these figures, one may witness the sensitiveness of buckling mode shapes towards scaling parameters. Also, these mode shapes help to predict the mechanical health and lifespan of several electromechanical devices.

Buckling Mode Shape
Buckling is the state of instability of structures that leads to structural failure. In this circumstance, the buckling mode shape has a vital role in predicting the instability. In this regard, buckling mode shapes were plotted with H a = 0.5 and L = 10 for different e 0 a (0.5, 1, 1.5, 2). Figures 7-10 show the buckling mode shapes for "P-P, C-P, C-C, and C-F" boundary conditions, respectively. From these figures, one may witness the sensitiveness of buckling mode shapes towards scaling parameters. Also, these mode shapes help to predict the mechanical health and lifespan of several electromechanical devices.

Buckling Mode Shape
Buckling is the state of instability of structures that leads to structural failure. In this circumstance, the buckling mode shape has a vital role in predicting the instability. In this regard, buckling mode shapes were plotted with for different a e 0 (0.5, 1, 1.5, 2). Figures 7-10 show the buckling mode shapes for "P-P, C-P, C-C, and C-F" boundary conditions, respectively. From these figures, one may witness the sensitiveness of buckling mode shapes towards scaling parameters. Also, these mode shapes help to predict the mechanical health and lifespan of several electromechanical devices.

Concluding Remarks
The buckling behavior of Electromagnetic nanobeam is investigated in the combined framework of Euler-Bernoulli beam theory and Eringen's nonlocal theory. Critical buckling load parameters were obtained using shifted Chebyshev polynomials-based Rayleigh-Ritz method for all the classical boundary conditions such as "Pined-Pined (P-P), Clamped-Pined (C-P), Clamped-Clamped (C-C), and Clamped-Free (C-F)". The C-F boundary condition converges faster with 5 = N , whereas other boundary conditions such as P-P, C-P, and C-C require 7 = N for achieving convergence to the desired accuracy. Critical buckling load parameters decrease with an increase in small scale parameter, and this decline is more in case of the C-C boundary condition. It may also be noted that the influence of small scale parameter is comparatively more in C-C edge and less in C-F edge. It is interesting to note that the critical buckling load decreases with an increase in aspect ratio. This decrease is more consequential for the lower values of aspect ratio. We may note that the critical buckling load decreases with an increase in Hartmann parameter, but this drop in critical buckling load is prolonged. The C-C nanobeam possesses the highest critical buckling load, whereas the C-F nanobeam possesses the lowest.