Dynamic stiffness formulation for a micro beam using Timoshenko–Ehrenfest and modified couple stress theories with applications

Several models within the framework of continuum mechanics have been proposed over the years to solve the free vibration problem of micro beams. Foremost amongst these are those based on non-local elasticity, classical couple stress, gradient elasticity and modified couple stress theories. Many of these models retain the basic features of the Bernoulli–Euler or Timoshenko–Ehrenfest theories, but they introduce one or more material scale length parameters to tackle the problem. The work described in this paper deals with the free vibration problems of micro beams based on the dynamic stiffness method, through the implementation of the modified couple stress theory in conjunction with the Timoshenko–Ehrenfest theory. The main advantage of the modified couple stress theory is that unlike other models, it uses only one material length scale parameter to account for the smallness of the structure. The current research is accomplished first by solving the governing differential equations of motion of a Timoshenko–Ehrenfest micro beam in free vibration in closed analytical form. The dynamic stiffness matrix of the beam is then formulated by relating the amplitudes of the forces to those of the corresponding displacements at the ends of the beam. The theory is applied using the Wittrick–Williams algorithm as solution technique to investigate the free vibration characteristics of Timoshenko–Ehrenfest micro beams. Natural frequencies and mode shapes of several examples are presented, and the effects of the length scale parameter on the free vibration characteristics of Timoshenko–Ehrenfest micro beams are demonstrated.


Introduction
Within the confines of continuum mechanics, there are many applications of micro beams, particularly in the design and developments of micro-electromechanical systems (Attia et al., 1998;Li et al., 2003;Moeenfard and Ahmadian, 2013;Pei et al., 2004). For such applications, classical beam theories are inapplicable because these theories ignore the microstructure dependent size effect. There are numerous experimental observations (Fleck et al., 1994;Chong and Lam, 1999;Chong et al., 2001;Lam et al., 2003;Ma and Clarke, 1995;Stolken and Evans, 1998) which strongly support that the size effect in a microstructure matters and its influence can be significant. The historical development of micro continuum mechanics took place more than a century ago, and the current state of the art owes a huge debt to early pioneers such as Cosserat and Cosserat (1909), Eringen (1966Eringen ( , 1972Eringen ( , 1983, Koiter (1964), Mindlin (1963Mindlin ( , 1964, Mindlin and Tiersten (1962) and Toupin (1962). Foremost amongst the micro continuum mechanical theories which survived the test of time and continues to receive attention even to this day are the non-local elasticity theory originally proposed by Eringen (1966Eringen ( , 1972Eringen ( , 1983 and later further developed by Reddy (2007) and Challamel and Wang (2008); the classical couple stress theory introduced by Koiter (1964), Mindlin and Tiersten (1962), Toupin (1962) and further used by Anthoine (2000); the gradient elasticity theory in connection with non-local elasticity presented by Mindlin (1963Mindlin ( , 1964 and subsequently applied by Akgoz and Civalek (2012), Askes and Aifantis (2011), Kong et al. (2009) and Lam et al. (2003);and the modified couple stress theory (MCST) developed by Yang et al. (2002) and further exploited by Attia and Emam (2018) and Park and Gao (2006). A survey paper by Thai et al. (2017) elucidates the developments of micro mechanics continuum models, and the paper provides considerable insights into the topic. Based on these micro continuum theories, the static and dynamic behaviour of micro beams with varying degrees of complexities has been investigated by many researchers in recent years (Ansari et al., 2015;Challamel et al., 2015;Darijini and Mohammadabadi, 2014;Essen, 2020;Fathalilou et al., 2014;Kong et al., 2008;Liang et al., 2014;Ma et al., 2008;Mustapha, 2020;Noori et al., 2016;Reddy, 2011;Yayli, 2017). For instance, within the chosen background theory of Bernoulli-Euler beam, Fathalilou et al. (2014) and Kong et al. (2008) used couple stress theory, whereas Challamel et al. (2015), Liang et al. (2014) and Yayli (2017) used strain-gradient based non-local elasticity theory to investigate the free vibration behaviour of micro beams. The subsequent developments for the free vibration analysis of micro beams followed quite intensely when the Timoshenko-Ehrenfest theory (Elishakoff, 2020a(Elishakoff, , 2020bTimoshenko, 1921Timoshenko, , 1922, which extends the Bernoulli-Euler beam theory by accounting for the effects of shear deformation and rotatory inertia, was implemented in the micro beam theory (Ansari et al., 2015;Ma et al., 2008;Reddy, 2011). It has come to the notice of the authors in recent years that Ehrenfest contributed in a substantial way to the enhancement of the Bernoulli-Euler beam theory to include the effects of shear deformation and rotatory inertia. In fact, the frequency equation of the so-called Timoshenko beam was jointly developed by Timoshenko and Ehrenfest before Timoshenko published his classic papers (Timoshenko, 1921(Timoshenko, , 1922, and the origin of the theory was recently researched by Elishakoff (2020aElishakoff ( , 2020b. The contribution of Ehrenfest was acknowledged by Timoshenko, see his footnote in Timoshenko (1921). Noori et al. (2016) published an interesting paper in which they developed higher order micro beam model to carry out the free vibration analysis. However, one of their conclusions 'The natural frequencies of beams with free and clamped edges employing Bernoulli-Euler model are identical, which is doubtful physically' is frankly debatable. Noori et al. (2016) appear to have overlooked that the beam with free-free boundary condition will have rigid-body modes of zero frequencies which might have been discounted by them when drawing this conclusion. Natural frequencies corresponding to the elastic modes for the free-free and clamped-clamed cases are equal, but such a comparison needs qualification based on the order of the natural frequencies.
In the current work, the Timoshenko-Ehrenfest theory in conjunction with the MCST which has the advantage of using only one material length scale parameter (Ma et al., 2008;Mustapha, 2020;Reddy, 2011;Yang et al., 2002) is a fundamental consideration to develop the dynamic stiffness theory and then carry out the free vibration analysis of the Timoshenko-Ehrenfest micro beam. In this context, it should be noted that linear small deflection theory has been used and any consideration of non-linear effects is outside the scope of this paper but interested readers are referred to Asghari et al., 2010;Farokhi et al., 2013;Ghayesh et al., 2013;Simsek, 2014;Wang et al., 2013 for further research on the subject.
Thus, the main purpose of this paper is to develop for the first time the dynamic stiffness matrix of a micro beam using Timoshenko-Ehrenfest theory and MCST and then carry out its free vibration analysis. The general procedure for dynamic stiffness theory development and its application can be found in the work of Banerjee (1997Banerjee ( , 2015, amongst others. In the current work, bending or flexural vibration is given precedence in the development of the dynamic stiffness theory because axial vibration is not affected by the small length parameter when the MCST is used (Ma et al., 2008;Mustapha, 2020;Reddy, 2011). The investigation is carried out in following steps. First, the potential and kinetic energies of a micro beam are formulated using Timoshenko-Ehrenfest theory and MCST. The former accounts for the shear deformation and rotatory inertia, whereas the latter accounts for the material scale length parameter. Thus, the formulation starts with the fundamental premises of MCST to obtain the expressions for the kinetic and potential energies, and then by using Hamilton's principle, the governing differential equations of motion in free vibration and natural boundary conditions of a micro beam based on Timosheno-Ehrenfest theory and MCST are derived. For harmonic oscillation, the governing differential equations are solved in closed analytical form. Next, the dynamic stiffness matrix of the micro beam based on the MCST and Timoshenko-Ehrenfest hypothesis is formulated by relating the amplitudes of forces and moments to those of the displacements and rotation at the ends of the micro beam element. Finally, the well-established Wittrick-Williams algorithm (Williams, 1993;Williams and Wittrick, 1983;Wittrick and Williams, 1971) is applied to the resulting dynamic stiffness matrix to compute the natural frequencies and mode shapes of some selected examples of micro beams. The effects of the small length scale parameter as well as those of shear deformation and rotatory inertia on the results are demonstrated, and some results are validated using published results. It is significant to note that the micro beam theory using MCST in combination with the Timoshenko-Ehrenfest theory gives rise to three displacement variables (excluding the axial displacement (u)) at each node of the beam in bending or flexure. These are bending or flexural displacement (v), bending or flexural rotation (θ) and the first derivative of the bending or flexural displacement (v 0 ). Unfortunately, the published literature on the free vibration of Timoshenk-Ehrenfest micro beams does not give due recognition to the displacement component v 0 in any practical or introspective manner. A secondary purpose of this paper is to redress this disparity.

Fundamental preliminaries
The MCST used here in the derivation of the governing differential equations of motion in free vibration of a micro beam in conjunction with the Timoshenko-Ehrenfest model was developed by Yang et al. (2002) with the underlying postulate that the strain energy is a function of both the strain tensor, which is conjugated to stress tensor, and importantly, the curvature tensor which is conjugated to couple stress tensor. Based on this postulate, the strain or potential energy U of a linearly deformed elastic body can be written in tensorial notation as (Kong et al., 2008;Park and Gao, 2006;Reddy, 2011) where σ ij and ε ij (i and j are the usually adopted dummy notations with i, j =1, 2, 3) are the components of the stress (σ) and strain (ε) tensors and m ij and χ ij signify the components of the deviatoric part of the symmetric couple stress tensor (m) and symmetric curvature tensor ðχÞ, respectively, and the integration is carried out over the entire volume of the elastic body. The expressions for σ ij , ε ij , m ij and χ ij in the usual notations are given by (Kong et al., 2008;Ma et al., 2008) where λ and G (often denoted by λ and μ) are Lame's constant, I is the unit or identity matrix (so that its elements δ ij is 1 when i = j and 0 when i ≠ j), u i and u j are the components of the displacement vector u, θ i and θ j are the components of rotation vector θ, and l is the material length scale parameter which identifies the effect of the couple stress as described by Yang et al. (2002). It should be noted that Chong et al. (2001) and Li et al. (2018) successfully determined the small length scale parameter l by experiment. For homogeneous and isotropic material, Lame's constants λ and G appearing in equations (2) and (4) can be found in standard texts in elasticity. These are where E is Young's modulus, ν is Poisson's ratio and G is the shear modulus or modulus of rigidity of the material.

Derivation of the governing differential equations
In a rectangular right-handed Cartesian coordinate system, Figure 1 shows a uniform micro beam of length L in a greatly exaggerated scale with the X-axis coinciding with the centroidal axis of the beam. The Y-axis is upward and perpendicular to the X-axis representing the direction of bending or flexural displacement of the beam. The Z-axis, which is not shown, is perpendicular to the plane of the paper represents the neutral axis of the beam. Based on the preliminaries of MCST (Ma et al., 2008;Mustapha, 2020;Reddy, 2011;Yang et al., 2002) briefly described in section 1.1, the strain or potential energy U of a linearly deformed elastic body like the micro beam shown in Figure 1 can now be formed, whereas the kinetic energy can be routinely expressed using the transverse and rotational components of a point lying on the neutral axis of the beam. Using the kinetic and potential energies formulated in this way, Hamilton's principle is applied here to derive the governing differential equations of motion of a Timoshenko-Ehrenfest micro beam undergoing free vibration.
Referring to Figure 1, let v be the bending or flexural displacement and θ be the bending rotation of the Timoshenko-Ehrenfest micro beam. The generalised displacement field of equations (3) and (5) then reduces to much simplified forms, as follows.
The bending rotation θðx,tÞ in equation (7) above is essentially the angle of rotation of the cross-section of the beam about the Z-axis which according to Timoshenko-Ehrenfest beam theory is given by   (8) where γ xy is the shearing strain (ε 12 in tensorial notation) which in fact is the slope arising from the transverse shear force that must be deducted from the total slope ð∂vðx,tÞ=∂xÞ of the flexural axis to give the bending slope or the angle of rotation θðx,tÞ of the cross-section. Equation (3) with the help of equations (7) and (8) gives the expression for the strains as The components of rotations based to the displacement field defined in equation (7) are given by With the help of equations (5) and (10), the components of the curvature tensor χ can now be obtained as The components of classical stress tensor σ and the deviatoric part of the couple stress tensor m from equations (2) and (4) are, respectively, given by Assuming linear small deflection theory, the expression for the potential energy U of the micro beam can be obtained by substituting equations (9)-(13) into equation (1) to give where A is the area, I is the second moment of area and k is the shear correction factor of the micro beam cross-section.
The expression of kinetic energy (T) is given by where ρ is density of the micro beam material and t is time. Substitution or equation (7) into Equation (15) gives The expressions for kinetic (T) and potential (U) energies developed above can now be used to derive the governing differential equations of motion of a micro beam undergoing free natural vibration by applying Hamilton's principle.
Hamilton's principle states where t 1 and t 2 are the time interval in the dynamic trajectory, and δ is the usual variational operator. Substituting T and U from equations (16) and (14) into equation (17) and then integrating by parts and collecting terms give the governing differential equations and natural boundary conditions as follows.
Governing differential equations kAG Natural boundary conditions Axial force Shear force Bending Moment Higher order moment Clearly, the axial motion is uncoupled from the bending or flexural motion, but importantly, the free vibration behaviour in axial motion is not affected by the material scale length parameter (l), unlike the case with bending of flexural vibration, see equations (18)-(20). Therefore, attention here is focused on bending or flexural vibration only when developing the dynamic stiffness method because the axial vibration being unaffected by material scale length parameter is covered by classical theory that can be found in stand texts.
However, it should be recognised that for a Timoshenko-Ehrenfest micro beam in bending or flexural motion, there are three displacement components v, θ and v 0 that arise will lead to a sixth order system, whereas in classical Timoshenko-Ehrenfest beam theory, the displacement field consists of v and θ only (without the v 0 ) resulting in a fourth order system. This difference was observed by earlier investigators (Ma et al., 2008;Mustapha, 2020;Reddy, 2011), but apparently, these investigators did not pay much attention to the additional displacement parameter v 0 , particularly when obtaining the results. One of the purposes of this investigation is to consider the new displacement component v 0 both in the theory as well as in the results.

Formulation of the dynamic stiffness matrix
We assume harmonic oscillation with circular or angular frequency ω rad/s for the Timoshenko-Ehrenfest micro beam undergoing free bending or flexural motion so that where V and Θ are the amplitudes of the bending displacement and bending rotation, respectively, and i ¼ ffiffiffiffiffiffi ffi À1 p .
Substituting equation (25) into equations (19) and (20) and introducing the non-dimensional length ξ = x/L led to the following ordinary differential equations It is now possible to eliminate either V or Θ from equations (26) and (27) to arrive at the following sixth order ordinary differential equation which is identically satisfied by both V and Θ. À and It should be noted that the differential equation for free vibration of a conventional beam of normal length using classical theory whether Bernoulli-Euler or Timoshenko-Ehrenfest theory results in a fourth order differential equation as opposed to the sixth order differential equation for a Timoshenko-Ehrenfest micro beam given by equation (28).
By assuming the solution in the form H ¼ e pξ where p is a constant, yet to be determined, the characteristic or auxiliary equation of the differential equation, equation (28) can be expressed as a 1 p 6 þ a 2 p 4 þ a 3 p 2 þ a 4 ¼ 0 The polynomial equation, equation (32) can be reduced to a cubic and then solved analytically using standard procedure (Pipes and Harvill, 2014). By taking the square root of each of the three roots of the cubic, which could be either real or complex, the six roots r j ðj ¼ 1; 2,…,6Þ of the characteristic or auxiliary equation (32) can be computed leading to the solutions of the differential equation, equation (28) as where R j and Q j ðj ¼ 1; 2,…,6Þ are two sets of different constants which can be related to each other, by using equations (26) and (27). The relationship between R j and Q j is obtained as where The constant vectors Q and R can be written as where the upper suffix T denotes a transpose.
The amplitude of the bending or flexural rotation in terms of R j from equations (33) and (34) is The first derivative of the amplitude of the bending or flexural displacement is The amplitudes of the shear force (S), bending moment (M) and higher order moment ðM Þ are obtained in terms of R j using equations (22)-(25) and equations (31), (33), (37) and (38) and after making some algebraic manipulation to give The dynamic stiffness matrix of the Timoshenko-Ehrenfest micro beam can now be formulated by applying the boundary conditions in algebraic form for both displacements and forces at the ends (nodes) of the beam. The procedure essentially eliminates the constants R j (j =1, 2 ,3….6) from equations (33) and (39) to construct the required dynamic stiffness matrix, providing the forcedisplacement relationship at the nodes. Referring to the sign convention for positive shear force, bending moment and higher order moment shown in Figure 2, the boundary conditions for displacements and forces are illustrated in Figure 3. These are: At Substituting equations (41) and (42) into equations (33), (37) and (38) gives the following matrix equation relating the displacements of the two ends with the constants R j (j =1, 2, 3, .. 6).
By eliminating the constant vector R from equations (44) and (46), F and δ can now be related to give the dynamic stiffness matrix relationship of the Timoshenko-Ehrenfest micro beam as is the required 6×6 frequency-dependent dynamic stiffness matrix. It should be noted that the resulting dynamic stiffness matrix K of equation (49) will be always symmetric and real as observed in an earlier investigation (Banerjee, 2003), with imaginary part of each element being zero although the matrices B and D are complex and asymmetric. The real part of the expanded dynamic stiffness matrix giving the relationship between the amplitudes of the forces to those of the displacements can be expressed in the following way The above frequency-dependent dynamic stiffness matrix K can now be used to compute the natural frequencies and mode shapes of either an individual Timoshenko-Ehrenfest micro beam or an assembly of them for different boundary conditions. A reliable and accurate method of solving the eigenvalue problem is to apply the Wittrick-Williams algorithm (Wittrick and Williams, 1971) which is well suited for the dynamic stiffness method applications. The algorithm uses the Sturm sequence property of the dynamic stiffness matrix and ensures that no natural frequencies of the structure analysed are missed. The detailed working procedure of the algorithm is not given here because it has featured in literally hundreds of papers in the literature. However, interested readers are referred to the original work of Wittrick and Williams (Wittrick and Williams, 1971) for a detailed insight.

Numerical results and discussion
The dynamic stiffness matrix developed above is applicable to free vibration analysis of a Timoshenko-Ehrenfest micro beam in bending or flexural motion, but the dynamic stiffness matrix in axial vibration (Howson et al., 1983) which is uncoupled from bending or flexural vibration can also be incorporated. However, it is well known that the small length parameter (l) does not affect the axial vibration characteristics of a micro beam when using the MCST (Ma et al., 2008;Mustapha, 2020;Reddy, 2011). Therefore, leaving aside the axial vibration, only the bending or flexural vibration is considered here when computing the numerical results.
The first set of results was computed for an individual micro beam of length L with rectangular cross-section of  width b and depth (or thickness) h, for three different boundary conditions, namely, the Simply supported-Simply supported (S-S), Clamped-Free (C-F) and Clamped-Clamped (C-C). The data used in the analysis are taken from (Reddy, 2011;Mustapha, 2020) and are given as follows.
The first five non-dimensional natural frequencies λ i ¼ ω i L 2 ffiffiffiffiffiffiffiffiffiffiffiffiffi ρA=EI p (i = 1, 2, ... 5) of the micro beam for the S-S, C-F and C-C boundary conditions are respectively shown in Tables 1-3 for various values of the size-dependent nondimensional material length scale parameter l/h (or l/L). Note that the value of l must not be literally zero in the data, but a negligibly small value close to zero, for example, l = h × 10 À6 can be used to avoid numerical overflow. From a computational accuracy standpoint, this is all right, bearing in mind that l = 0 corresponds to the case of classical Bernoulli-Euler or Timoshenko-Ehrenfest beam for which results can be found in standard texts.
Comparative results for the first three natural frequencies of the Timoshenko-Ehrenfest micro beam for S-S boundary condition have been published by Reddy (2011) and are shown in Table 1 alongside the results computed using the present theory. Clearly, the results computed using the present theory are in excellent agreement with those of Reddy (2011).
The natural frequencies for C-F boundary condition of the micro beam are shown in Table 2. The first three natural frequencies for this problem for l/h = 0.0 and 1.0 published recently by Mustapha (2020) who used finite element method, but importantly, accounted for the shear deformation and rotatory inertia effects as in the present case, agreed very well with the above results, but are not shown for brevity, the disagreement being less than 5%.
For completeness, the natural frequencies of the Timoshenko-Ehrenfest micro beam for the C-C boundary condition were also computed, and they are shown in Table 3. The authors could not find suitable comparable results in the literature for this case.
Representative mode shapes of the Timoshenko-Ehrenfest micro beam for S-S, C-F and C-C boundary conditions with l=h ¼ 0:6 are illustrated in Figures 4-6, respectively, showing bending displacement (V), bending rotation (Θ) and the first derivative of the bending displacement ðV 0 Þ. To be consistent with the units, the bending rotation (Θ) and the first derivative of the bending displacement ðV 0 Þ were multiplied by the width (b) of the Timoshenko-Ehrenfest micro beam when plotting the mode shapes. It is clear from the mode shapes that for lower order modes, the bending rotation (Θ) and the first derivative of the bending displacement ðV 0 Þ are almost indistinguishable, but for higher order modes, some differences in their values exist.
The final set of results was computed for a stepped micro beam with cantilever boundary condition shown in Figure 7 by using the above dynamic stiffness theory. Each of the two segments of the stepped micro beam has a solid circular cross-section with diameters d 1 = 25 × 10 À6 m and d 2 = 15 × 10 À6 m, respectively, with the corresponding lengths of the two segments are taken as L 1 = 100 × 10 À6 m and L 2 = 125 × 10 À6 m so that the total length L = 225 × 10 À6 m. Young's modulus (E), density (ρ), Poisson's ratio (ν) of material and shear correction factor k are taken to be the same as the previous examples. Table 4 shows the first five natural frequencies of the stepped micro beam for different values of l=L (or l=d 2 ) ratios using both Timoshenko-Ehrenfest and Bernoulli-Euler dynamic theories. It should be noted that the results for l=L ¼ 0 correspond to classical Bernoulli-Euler and/or Timoshenko-Ehrenfest theories which do account for the small-scale effect. These results for l=L ¼ 0 when using Bernoulli-Euler and Timoshenko-Ehrenfest dynamic theories were further checked using some well-established computer programs (Anderson et al., 1986;Howson et al., 1983). Clearly, for higher values of l=L (or l=d 2 ) and for higher natural frequencies, the differences in the results using the two theories are pronounced, but importantly, the effects of shear deformation and rotatory inertia are more pronounced in micro beams when compared to classical beams. The Timoshenko-Ehrenfest theory gives lower values than that of the Bernoulli-Euler theory as expected. The percentage difference between the two theories in the third natural frequency for different values of l=L is shown in Figure 8.

Conclusions
Using the Timoshenko-Ehrenfest beam theory and the modified couple stress model, the dynamic stiffness matrix of a micro beam is developed. The ensuing dynamic stiffness matrix which relates the amplitudes of the forces to  those of the corresponding displacements of the harmonically vibrating Timoshenko-Ehrenfest micro beam is applied in conjunction with the Wittrick-Williams algorithm to compute the natural frequencies and mode shapes for different boundary conditions of some illustrative examples. A substantial part of the results is validated against published results. To demonstrate the versatility of the method, a stepped micro beam is also analysed. One of the displacement fields which is the first derivative of the flexural displacement appearing in the governing differential equations and natural boundary conditions has not been given due recognition in the literature. The current investigation has addressed this issue and provided the influence of this displacement on results. The results given in the paper can be used as benchmark to validate finite element and other approximate methods.