Variational Formulations and Isogeometric Analysis of Timoshenko–Ehrenfest Microbeam Using a Reformulated Strain Gradient Elasticity Theory

Abstract

This paper presents a novel non-classical Timoshenko–Ehrenfest beam model based on a reformulated strain gradient elasticity theory. The strain gradient effect, couple stress effect, and velocity gradient effect for vibration are included in the new model by only one material length scale parameter for each. The variational formulation and Hamilton’s principle are applied to derive the governing equations and boundary conditions. Both an analytical solution and an isogeometric analysis approach are proposed for static bending and free vibration of the microbeam. A non-uniform rational B-splines (NURBS) isogeometric analysis with high-order continuity can effectively fulfill the higher derivatives of the displacement variables in the reformulated gradient beam model. Convergence studies and comparisons to the corresponding analytical solutions verify the model’s performance and accuracy. Finally, different boundary conditions, material length scale parameters, and beam thicknesses are investigated in order to certify the applicability of the proposed approach.

1. Introduction

Micro- and nanobeam structures are widely used in micro/nano-scale devices and MEMS such as atomic force microscopes, sensors, and actuators. However, as microstructure-dependent size effects [1,2] exist in these micro/nanobeams, classical continuum mechanics models without any material length scale parameter fail to describe the size effects [3,4]. To overcome this deficiency, several higher-order continuum theories have been developed to capture size effects, including strain gradient elasticity theory [5], nonlocal elasticity theory [6,7], and couple stress elasticity theory [8,9,10].

In a higher-order continuum theory, microstructure effects are considered via microstructure-dependent material parameters. The strain gradient elasticity theory (or Eringen–Mindlin micromorphic theory) [5,11,12,13,14], for instance, contains sixteen microstructure-dependent material parameters. This general theory appears in three forms with five additional material scale parameters for isotropic and centrosymmetric materials in the strain energy density function [15,16]. As the large number of material parameters are difficult to obtain, several simpler variants of strain gradient elasticity theory have been developed. Inspired by Fleck and Hutchinson [17], Lam et al. [1] presented a modified strain gradient elasticity theory (MSGET) reduced to three additional material parameters for size effects. Sedighi et al. [18] utilized strain gradient elasticity theory to analyze the size-dependent electromechanical instability of a cantilever nanoactuator. Kong et al. [19] utilized a modified strain gradient elasticity theory to study Euler–Bernoulli microbeams, while Abbasi [20] applied it to study size-dependent vibration problems in an atomic force microscope with an assembled cantilever probe (ACP). Wang et al. [21] studied the bending and vibration behavior of porous metal foam microbeams based on sinusoidal beam theory and MSGET. Based on Form II of Mindlin’s strain gradient theory [15,22], Altan and Aifantis [23] proposed a simplified strain gradient elasticity theory (SSGET) [24,25]. This simplified strain gradient elasticity theory includes only one additional material parameter, and is known as the dipolar gradient elasticity theory [26]. The SSGET has been utilized to investigated a variety of microstructure beam problems [27,28]. Khakalo et al. [29] researched the static bending, buckling and vibration behavior of 2D triangular lattices utilizing SSGET. Ansari and Torabi [30] studied the vibration behavior of carbon nanocones (CNCs) embedded in an elastic foundation using the variational differential quadrature method, which they actualized based on nonlocal elasticity theory [6,7]. On the basis of the same theory, Ansari et al. [31] studied vibration problems in circular double-layered graphene sheets (DLGSs) under thermal load and elastic foundation. Based on nonlocal strain gradient theory, Torabi and Zabihi et al. [32,33] developed an analytical nanoplate model for the static and dynamic pull-in instability of functionally graded nanoplates. On the other hand, the classical couple stress elasticity theory [5,8,10] was developed using an additional two material parameters, while the modified couple stress theory (MCST) [34,35] with only one microstructure-dependent material parameter was developed by ignoring the symmetric second gradient of displacement. Subsequently, a great number of microstructure beam models based on MCST have been presented [36,37,38,39,40,41,42,43].

However, the above works cannot consider the strain gradient and couple stress effects simultaneously. For vibration analysis, Papargyri-Beskou et al. [25] studied the velocity gradient effect for wave propagation using simplified strain gradient elasticity theory. Recently, Zhang and Gao [44] proposed a reformulated strain gradient elasticity theory which simultaneously considered the couple stress, strain gradient, and velocity gradient effects for vibration. On the other hand, in this theory only one material length scale parameter is needed for each of the couple stress, strain gradient, and velocity gradient effect. In addition, it can be easily degenerated into the MCST. In this work, the reformulated strain gradient elasticity theory is extended to develop a non-classical Timoshenko–Ehrenfest beam model. Based on this Timoshenko–Ehrenfest model, variational formulation, and Hamilton’s principle, the static bending and free vibration of microbeams can be solved by analytical solution and isogeometric analysis, effectively fulfilling the higher derivatives of the displacement variables in strain gradient elasticity theory.

Isogeometric analysis (IGA) is a new higher-order numerical approach developed by Hughes et al. [45] which has been widely applied in many fields of engineering, such as beam, plate, and shell structure analysis [46,47,48,49,50], contact analysis [51], damage and fracture mechanics [52,53,54], electromagnetic analysis [55], fluid mechanics [56], fluid–structure interaction [57], structure optimization [58,59], and more. Because of the high-order continuity requirement of generalized displacements in strain gradient elasticity theory, FEM with addition of a number of unknowns [60,61] is applied to obtain the second derivative of displacements, which is computationally expensive and unstable [62,63]. As the IGA can effectively satisfy the high-order continuity requirements, it is more suitable for investigating the microstructure-dependent behavior of microbeam/plate structures based on non-local elasticity theory [64], modified strain gradient elasticity theory [65,66], and MCST [67,68,69]. For example, Niiranen et al. [70,71,72] studied the fourth-order and sixth-order boundary value problems of bar and plane strain/stress models for static and vibration problem with IGA. Natarajan et al. [73] discussed the scale-dependent linear free flexural vibration behavior of functionally graded (FG) nanoplates using IGA-based gradient elasticity theory [74]. Based on the user element of Abaqus finite-element software, an isogeometric analysis of higher-order strain gradient elasticity was actualized by Khakalo and Niiranen [75]. Greco and Cuomo [76] presented an isogeometric approach for the numerical analysis of 3D Kirchhoff–Love rod theory. Balobanov et al. [77] developed an IGA for size-dependent analysis of a Timoshenko–Ehrenfest beam utilizing strain gradient elasticity theory. Based on the same theory, Tran and Niiranen [78] founded a nonlinear Euler–Bernoulli beam model applied to a lattice structure with IGA. Niiranen et al. [79] proposed an IGA for Euler–Bernoulli microbeams. Yaghoubi et al. [80] utilized IGA for the dynamics of anisotropic Bernoulli–Euler shear-deformable beams. However, all of the above works consider only the couple stress or strain gradient effect.

This work aims to develop a non-classical Timoshenko–Ehrenfest beam model based on a reformulated strain gradient elasticity theory. The strain gradient, couple stress, and velocity gradient effects (for vibration) are included by only one material parameter for each. Both an analytical solution and an isogeometric analysis approach are proposed in order to study the size effects of the Timoshenko–Ehrenfest microbeam model. The rest of this paper is structured as follows. Section 2 briefly introduces the theoretical formulation of the reformulated strain gradient elasticity theory. The variational formulation, motion equations, and boundary conditions for the non-classical Timoshenko–Ehrenfest beam model are described in Section 3. Section 4 focuses on the discretization equation of the NURBS-based isogeometric analysis. In Section 5, the static bending and free vibration of the Timoshenko–Ehrenfest microbeam are studied and the results obtained with the analytical and isogeometric analysis solutions are compared. Finally, the conclusions are provided in Section 6.

2. Reformulated Strain Gradient Elasticity Theory

This section provides a brief introduction to the reformulated strain gradient elasticity theory (RSGT) proposed by Zhang and Gao [44]. The RSGT incorporates two additional material parameters to describe the strain gradient and couple stress effects, respectively. When referring to RSGT, the total strain energy density function w is written as

$$w({\epsilon }_{ij},{\eta }_{ijk}^s,{\chi }_{ij})=\frac{1}{2}\lambda {\epsilon }_{ii}{\epsilon }_{jj}+\mu {\epsilon }_{ij}{\epsilon }_{ij}+l_s^2\mu {\eta }_{ijk}^s{\eta }_{ijk}^s+\mu l_m^2{\chi }_{ij}{\chi }_{ij}$$

(1)

where l_s and l_m are material scale parameters depicting the strain gradient effect and the couple stress effect, respectively; ${\epsilon }_{ij},\hspace{0.17em}{\chi }_{ij}$ and ${\eta }_{ijk}^s$ are the strain tensor, the curvature tensor, and the symmetric part of the second-order displacement gradient tensor, respectively. Two classical Lamé constants, λ and μ, are obtained as

$$\lambda =\frac{Ev}{(1+v)(1-2v)}$$

(2)

$$\mu =\frac{E}{2(1+v)}$$

(3)

where E and v are the Young’s modulus and Poisson ratio, respectively.

Based on Equation (1), the isotropic linear elastic material total strain energy U in region Ω is obtained as

$$U={\displaystyle {\int }_{\Omega }wdV=\frac{1}{2}}{\displaystyle {\int }_{\Omega }({\sigma }_{ij}{\epsilon }_{ij}+{\tau }_{ijk}^s{\eta }_{ijk}^s+m_{ij}{\chi }_{ij})}dV$$

(4)

where ${\sigma }_{ij}$, $m_{ij}$, and ${\tau }_{ijk}^s$ are respectively the Cauchy stress, the couple stress, and the symmetric part of the double stress tensors obtained by

$${\sigma }_{ij}=\lambda {\epsilon }_{kk}{\delta }_{ij}+2\mu {\epsilon }_{ij}$$

(5)

$${\tau }_{ijk}^s=2l_s^2\mu {\eta }_{ijk}^s$$

(6)

$$m_{ij}=2l_m^2\mu {\chi }_{ij}$$

(7)

where ${\delta }_{ij}$ refers to the Kronecker delta; ${\epsilon }_{ij}$, ${\chi }_{ij}$, and ${\eta }_{ijk}^s$ can be expressed as

$${\epsilon }_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})$$

(8)

$${\eta }_{ijk}^s=\frac{1}{3}(u_{i,jk}+u_{j,ki}+u_{k,ij})$$

(9)

$${\chi }_{ij}=\frac{1}{2}({\theta }_{i.j}+{\theta }_{j,i})$$

(10)

where $u_i$(i∈{x, y, z}) is the displacement vector and ${\theta }_i$ is the rotation vector, defined by

$${\theta }_i=\frac{1}{2}{e}_{ijk}{u}_{k,j}$$

(11)

where $e_{ijk}$ denotes the Levi–Civita symbol.

3. Non-Classical Timoshenko–Ehrenfest Beam Model

Consider a Timoshenko–Ehrenfest beam with a uniform cross section, as shown in Figure 1; the beam is subjected to a static load q(x) distributed along the x axis within the Cartesian coordinate system (x, y, z). The displacement field of Timoshenko–Ehrenfest beam theory can be depicted as [81,82]

$$u_1=-z\phi (x,t)$$

(12)

$$u_2=0$$

(13)

$$u_3=w(x,t)$$

(14)

where $u_i$(i = 1, 2, 3) denote the displacement components of a discretional point (x, y, z) on the cross-section beam, w represents the z-component of the displacement on the centroidal axis (i.e., the x-axis), and φ represents the angle of rotation about the y-axis.

By substituting Equations (12)–(14) into Equation (8), the strain tensor, ${\epsilon }_{ij}$, is obtained as

$${\epsilon }_{xx}=-z\frac{\partial \phi }{\partial x}, {\epsilon }_{xz}=\frac{1}{2}\left(\frac{\partial w}{\partial x}-\phi \right), {\epsilon }_{yy}={\epsilon }_{zz}={\epsilon }_{xy}={\epsilon }_{yz}=0$$

(15)

From Equations (11) and (12)–(14) we can find that

$${\theta }_y=-\frac{1}{2}\left(\frac{\partial w}{\partial x}+\phi \right), {\theta }_x={\theta }_z=0$$

(16)

By inserting Equation (16) into Equation (10), the components of the curvature tensor are obtained as

$${\chi }_{xy}=-\frac{1}{4}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right), {\chi }_{xx}={\chi }_{yy}={\chi }_{zz}={\chi }_{xz}={\chi }_{yz}=0$$

(17)

Substituting Equations (12)–(14) into Equation (9), the components of the second-order displacement gradient tensor can be written as

$$\begin{gathered} {\eta }_{xxx}^s=-z\frac{{\partial }^2\phi }{\partial x^2}, {\eta }_{xxz}^s=\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}-\frac{2}{3}\frac{\partial \phi }{\partial x}, \\ {\eta }_{yyy}^s={\eta }_{zzz}^s={\eta }_{xxy}^s={\eta }_{xyy}^s={\eta }_{xyz}^s={\eta }_{xzz}^s={\eta }_{yyz}^s={\eta }_{yzz}^s=0 \end{gathered}$$

(18)

Note that ${\eta }_{xxx}^s$ is of a higher order than the second-order micro displacement gradient, ${\eta }_{xxz}^s$ [83]. Hence, in the rest of the current formulation, ${\eta }_{xxx}^s$ can be neglected.

Substituting Equations (2), (3) and (15) into Equation (5) yields [38]

$$\begin{gathered} {\sigma }_{xx}=-\frac{E\left(1-v\right)}{\left(1+v\right)\left(1-2v\right)}z\frac{\partial \phi }{\partial x}, {\sigma }_{yy}={\sigma }_{zz}=-\frac{Ev}{(1+v)(1-2v)}z\frac{\partial \phi }{\partial x}, \\ {\sigma }_{xz}=\frac{E}{2(1+v)}\left(\frac{\partial w}{\partial x}-\phi \right), {\sigma }_{xy}={\sigma }_{yz}=0 \end{gathered}$$

(19)

Using Equations (7) and (17), the components of the modified couple stress tensor, $m_{ij}$, are obtained as

$$m_{xy}=-\frac{1}{2}\mu l_m^2\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right), m_{xx}=m_{yy}=m_{zz}=m_{xz}=m_{yz}=0$$

(20)

By substituting Equation (18) into Equation (6), the components of the double stress tensor, ${\tau }_{ijk}^s$, can be described as

$${\tau }_{xxz}^s=2l_s^2\mu \left(\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}-\frac{2}{3}\frac{\partial \phi }{\partial x}\right), {\tau }_{yyy}^s={\tau }_{zzz}^s={\tau }_{xxy}^s={\tau }_{xyy}^s={\tau }_{xyz}^s={\tau }_{xzz}^s={\tau }_{yyz}^s={\tau }_{yzz}^s=0$$

(21)

Using Equation (4) and Equations (19)–(21), the first variation of the strain energy U over time span [0, T] is obtained as

$$\delta {\displaystyle {\int }_0^{T}Udt}={\displaystyle {\int }_0^T{\displaystyle {\int }_{\Omega }\left({\sigma }_{xx}\delta {\epsilon }_{xx}+2{\sigma }_{xz}\delta {\epsilon }_{xz}+2m_{xy}\delta {\chi }_{xy}+3{\tau }_{xxz}^s\delta {\eta }_{xxz}^s\right)}dVdt}$$

(22)

Substituting Equations (15)–(18) into Equation (22) provides

$$\begin{array}{l}\delta {\displaystyle {\int }_0^{T}Udt={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L[-M_x\delta \frac{\partial \phi }{\partial x}+R_x\delta \left(\frac{\partial w}{\partial x}-\phi \right)-\frac{1}{2}{Y}_x\delta \left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right)}}}+Q_x\delta \left(\frac{{\partial }^{2}w}{\partial x^2}-2\frac{\partial \phi }{\partial x}\right)]dxdt\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left(\frac{\partial M_x}{\partial x}-R_x+\frac{1}{2}\frac{\partial Y_x}{\partial x}+2\frac{\partial Q_x}{\partial x}\right)\delta \phi dxdt}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left(-\frac{\partial R_x}{\partial x}-\frac{1}{2}\frac{{\partial }^2{Y}_x}{\partial x^2}+\frac{{\partial }^2{Q}_x}{\partial x^2}\right)}}\delta wdxdt+\hspace{0.17em}{\displaystyle {\int }_0^T[(-\frac{1}{2}{Y}_x-M_x-2}{Q}_x)\delta \phi \hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}\left(R_x+\frac{1}{2}\frac{\partial Y_x}{\partial x}-\frac{\partial Q_x}{\partial x}\right)\delta w+\left(Q_x-\frac{1}{2}{Y}_x\right)\delta {\frac{\partial w}{\partial x}]|}_0^{L}dt\end{array}$$

(23)

where

$$M_x={\displaystyle {\int }_A{\sigma }_{xx}zdA,\hspace{0.17em}\hspace{0.17em}{R}_x={\displaystyle {\int }_A{\sigma }_{xz}dA,\hspace{0.17em}\hspace{0.17em}{Y}_x={\displaystyle {\int }_A{m}_{xy}dA}},}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{Q}_x={\displaystyle {\int }_A{\tau }_{xxz}^{s}dA}$$

(24)

In Equation (24), the stress resultants L and A are, respectively, the beam length and the beam cross-sectional area.

The kinetic energy of the current Timoshenko–Ehrenfest beam can be described as [5,25]

$$K={\displaystyle {\int }_{\Omega }\frac{1}{2}}\rho \left[\frac{\partial u_i}{\partial t}\frac{\partial u_i}{\partial t}+l_v^2\frac{{\partial }^2{u}_i}{\partial x_j\partial t}\frac{{\partial }^2{u}_i}{\partial x_j\partial t}\right]dV$$

(25)

where $\rho$ and $l_v$ are the material density and the velocity gradient coefficient, respectively.

Using Equations (12)–(14) and (25), the first variation of kinetic energy can be rewritten as

$$\begin{array}{l}\delta {\displaystyle {\int }_0^{T}Kdt={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L[m_2\frac{\partial \phi }{\partial t}\delta \frac{\partial \phi }{\partial t}+m_0\frac{\partial w}{\partial t}}}}\delta \frac{\partial w}{\partial t}+l_v^2{m}_2\frac{{\partial }^2\phi }{\partial x\partial t}\delta \frac{{\partial }^2\phi }{\partial x\partial t}\hspace{0.17em}+l_v^2{m}_0\frac{\partial \phi }{\partial t}\delta \frac{\partial \phi }{\partial t}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+l_v^2{m}_0\frac{{\partial }^{2}w}{\partial x\partial t}\delta \frac{{\partial }^{2}w}{\partial x\partial t}]dxdt\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=-{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L[\left(m_2\frac{{\partial }^2\phi }{\partial t^2}-l_v^2{m}_2\frac{{\partial }^4\phi }{\partial x^2\partial t^2}+m_0{l}_v^2\frac{{\partial }^2\phi }{\partial t^2}\right)\delta \phi +m_0\delta w\left(\frac{{\partial }^{2}w}{\partial t^2}-l_v^2\frac{{\partial }^{4}w}{\partial x^2\partial t^2}\right)]}}dxdt\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{{\displaystyle {\int }_0^T\left(l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}\delta \phi +l_v^2{m}_0\frac{{\partial }^{3}w}{\partial x\partial t^2}\delta w\right)}|}_0^{L}dt\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}$$

(26)

where

$$m_0={\displaystyle {\int }_A\rho dA}$$

(27)

$$m_2={\displaystyle {\int }_A\rho z^{2}dA}$$

(28)

It is noteworthy that ρ as a constant over the cross-section is independent of time t, which results in ${\dot{m}}_0=0$ and ${\dot{m}}_2=0$ over time interval [0, T].

Based on modified couple stress theory (MCST) [38] and strain gradient theory [5], the first variations of the work by external forces on the current Timoshenko–Ehrenfest beam can be described as

$$\delta {\displaystyle {\int }_0^{T}Wdt}={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L(q\delta w+c\delta {\theta }_y)dxdt}}+{{\displaystyle {\int }_0^T\left[\left(\overline{N}\delta w-\overline{M}\delta \phi -\overline{H}\delta \frac{\partial w}{\partial x}\right)\right]}|}_0^{L}dt$$

(29)

where q and c are the z-component of body force and the y-component of body couple per unit length along the x-axis, respectively, ${\theta }_y$ is provided by Equation (16), $\overline{N}$ and $\overline{M}$ are the transverse force and the bending moment, respectively, and $\overline{H}$ is the high-order bending moment.

Based on Hamilton’s principle [38,84],

$$\delta {\displaystyle {\int }_0^T\left[K-\left(U-W\right)\right]}dt=0$$

(30)

By substituting Equations (23), (26) and (29) into Equation, (29) we can obtain

$$\begin{array}{l}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L-[\left(m_0\frac{{\partial }^{2}w}{\partial t^2}-m_0{l}_v^2\frac{{\partial }^{4}w}{\partial x^2\partial t^2}-\frac{\partial R_x}{\partial x}-\frac{1}{2}\frac{{\partial }^2{Y}_x}{\partial x^2}+\frac{{\partial }^2{Q}_x}{\partial x^2}-q-\frac{1}{2}\frac{\partial c}{\partial x}\right)\delta w}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left(m_2\frac{{\partial }^2\phi }{\partial t^2}-m_2{l}_v^2\frac{{\partial }^4\phi }{\partial x^2\partial t^2}+m_0{l}_v^2\frac{{\partial }^2\phi }{\partial t^2}+\frac{\partial M_x}{\partial x}-R_x+\frac{1}{2}\frac{\partial Y_x}{\partial x}+2\frac{\partial Q_x}{\partial x}+\frac{1}{2}c\right)\delta \phi ]dxdt\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-{\displaystyle {\int }_0^T[\left(l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}+R_x+\frac{1}{2}\frac{\partial Y_x}{\partial x}-\frac{\partial Q_x}{\partial x}-\overline{N}+\frac{1}{2}c\right)\delta w+}(Q_x-\frac{1}{2}{Y}_x+\overline{H})\delta \frac{\partial w}{\partial x}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\left(l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}-M_x-\frac{1}{2}{Y}_x-2Q_x+\overline{M}\right)\delta \phi ]|}_0^{L}dt=0\end{array}$$

(31)

Considering the arbitrariness of $\delta \phi$ and $\delta w$ for a given x ∈ [0, L] and t ∈ [0, T] and applying the essential law of the calculus of variations [85] in Equation (31) results in the following equations of motion:

$$m_0\left(\frac{{\partial }^{2}w}{\partial t^2}-l_v^2\frac{{\partial }^{4}w}{\partial x^2\partial t^2}\right)-\frac{\partial R_x}{\partial x}-\frac{1}{2}\frac{{\partial }^2{Y}_x}{\partial x^2}+\frac{{\partial }^2{Q}_x}{\partial x^2}=q+\frac{1}{2}\frac{\partial c}{\partial x}$$

(32)

$$m_2\left(\frac{{\partial }^2\phi }{\partial t^2}-l_v^2\frac{{\partial }^{4}w}{\partial x^2\partial t^2}\right)+m_0{l}_v^2\frac{{\partial }^2\phi }{\partial t^2}+\frac{\partial M_x}{\partial x}-R_x+\frac{1}{2}\frac{\partial Y_x}{\partial x}+2\frac{\partial Q_x}{\partial x}+\frac{1}{2}c=0$$

(33)

with the complete boundary conditions

$$l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}+R_x+\frac{1}{2}\frac{\partial Y_x}{\partial x}-\frac{\partial Q_x}{\partial x}+\frac{1}{2}c=\overline{N} \mathrm{or} w=\overline{w} \mathrm{at} x=0 \mathrm{and} x=L$$

(34)

$$M_x+\frac{1}{2}{Y}_x+2Q_x-l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}=\overline{M} \mathrm{or} \phi =\overline{\phi } \mathrm{at} x=0 \mathrm{and} x=L$$

(35)

$$\frac{1}{2}{Y}_x-Q_x=\overline{H} \mathrm{or} \frac{\partial w}{\partial x}=\overline{\frac{\partial w}{\partial x}} \mathrm{at} x=0 \mathrm{and} x=L$$

(36)

where the overhead bar denotes the prescribed value.

Inserting Equations (19)–(21) into Equation (24) yields

$$M_x=-\frac{E(1-v)I}{(1+v)(1-2v)}\frac{\partial \phi }{\partial x}$$

(37)

$$R_x=K_S\mu A\left(\frac{\partial w}{\partial x}-\phi \right)$$

(38)

$$Y_x=-\frac{1}{2}\mu l_m^{2}A\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right)$$

(39)

$$Q_x=2l_s^2\mu A\left(\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}-\frac{2}{3}\frac{\partial \phi }{\partial x}\right)$$

(40)

where K_S is the shear coefficient, defined as Ks = (5 + 5v)/(6 + 5v) [38], and I is the moment of inertia of the cross-sectional area about the y-axis, obtained as

$$I={\displaystyle {\int }_A{z}^{2}dA}$$

(41)

Substituting Equations (34)–(40) into Equations (32) and (33), the equations of motion for Timoshenko–Ehrenfest beam with respect to $\phi$ and w can be described as

$$\begin{array}{l}{m}_0\left(\frac{{\partial }^{2}w}{\partial t^2}-l_v^2\frac{{\partial }^{4}w}{\partial x^2\partial t^2}\right)-K_s\mu A\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+\frac{1}{4}\mu l_m^{2}A\left(\frac{{\partial }^{4}w}{\partial x^4}+\frac{{\partial }^3\phi }{\partial x^3}\right)\\ \hspace{0.17em}\hspace{0.17em}=q+\frac{1}{2}c-2l_s^2\mu A\left(\frac{1}{3}\frac{{\partial }^{4}w}{\partial x^4}-\frac{2}{3}\frac{{\partial }^3\phi }{\partial x^3}\right)\end{array}$$

(42)

$$\begin{array}{l}\frac{E\left(1-v\right)I}{\left(1+v\right)\left(1-2v\right)}\frac{{\partial }^2\phi }{\partial x^2}+K_s\mu A\left(\frac{\partial w}{\partial x}-\phi \right)+\frac{1}{4}\mu l_m^{2}A\left(\frac{{\partial }^{3}w}{\partial x^3}+\frac{{\partial }^2\phi }{\partial x^2}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=m_2\left(\frac{{\partial }^2\phi }{\partial t^2}-l_v^2\frac{{\partial }^4\phi }{\partial x^2\partial t^2}\right)+m_0{l}_v^2\frac{{\partial }^2\phi }{\partial t^2}+\frac{1}{2}c+4l_s^2\mu A\left(\frac{1}{3}\frac{{\partial }^{3}w}{\partial x^3}-\frac{2}{3}\frac{{\partial }^2\phi }{\partial x^2}\right)\end{array}$$

(43)

It is clearly seen from Equations (42) and (43) that the present Timoshenko–Ehrenfest beam model incorporates three additional material parameters (i.e., $l_s$, $l_m$, and $l_v$), which enables the current model to describe the size-dependent elastic properties of the microstructure.

When $l_s$ = $l_v$ = 0, the governing equations from Equations (42) and (43) degenerate to the modified couple stress theory formulations [38], as follows:

$$m_0\frac{{\partial }^{2}w}{\partial t^2}-K_s\mu A\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+\frac{1}{4}\mu l_m^{2}A\left(\frac{{\partial }^{4}w}{\partial x^4}+\frac{{\partial }^3\phi }{\partial x^3}\right)=q+\frac{1}{2}\frac{\partial c}{\partial x}$$

(44)

$$\frac{E\left(1-v\right)I}{\left(1+v\right)\left(1-2v\right)}\frac{{\partial }^2\phi }{\partial x^2}+K_s\mu A\left(\frac{\partial w}{\partial x}-\phi \right)+\frac{1}{4}\mu l_m^{2}A\left(\frac{{\partial }^{3}w}{\partial x^3}+\frac{{\partial }^2\phi }{\partial x^2}\right)-\frac{1}{2}c=m_2\frac{{\partial }^2\phi }{\partial t^2}$$

(45)

When $l_s$ = $l_m$ = $l_v$ = 0 and c = 0, the governing equations from Equations (42) and (43) reduce to the classical Timoshenko–Ehrenfest beam model, which can be obtained as

$$m_0\frac{{\partial }^{2}w}{\partial t^2}=K_s\mu A\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+q$$

(46)

$$m_2\frac{{\partial }^2\phi }{\partial t^2}=EI\frac{{\partial }^2\phi }{x^2}+K_s\mu A\left(\frac{\partial w}{\partial x}-\phi \right)$$

(47)

Furthermore, the equilibrium equations and boundary conditions from Equations (34)–(36), (42) and (43), with no shear deformation, reduce to a micro-scale Bernoulli–Euler beam based on reformulated strain gradient elasticity theory [44].

4. Isogeometric Analysis Approach

This section is devoted to the development of an isogeometric analysis (IGA) approach based on the above non-classical Timoshenko–Ehrenfest beam model for static bending and free vibration analysis. The numerical results predicted via both analytical solution and isogeometric analysis are presented in the following section.

4.1. NURBS Basis Functions

For a straight beam, a NURBS curve can be degenerated into a B-spline curve and NURBS basis functions with a polynomial order p, which consist of weight B-spline functions. These can be described as

$$R_{i,p}(\xi )=\frac{N_{i,p}(\xi )w_i}{{\displaystyle \sum _{j=1}^n{N}_{j,p}(\xi )}{w}_j}$$

(48)

where $w_i$ and $\xi$ are the $i^{th}$ weight and the parametric coordinate, respectively, n depicts the number of control points and NURBS basis functions, and $N_{i,p}(\xi )$ is the $i^{th}$ B-splines basis function of degree p, obtained as

$$N_{i,0}(\xi )=\{\begin{array}{c}\hspace{0.17em}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{if}\hspace{0.17em}{\xi }_i\le \xi \le {\xi }_{i+1}\\ \hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{otherwise}\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p=0$$

(49)

and

$$N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}\left(\xi \right)+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}\left(\xi \right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{for}\hspace{0.17em}\hspace{0.17em}p\ge 1$$

(50)

4.2. Discretized Equations

Based on the current Timoshenko–Ehrenfest beam model, the NURBS basis functions are used to approximate the generalized mid-plane displacement, as follows:

$${\mathbf{u}}^h={\displaystyle \sum _I^{NP}{R}_I{\mathbf{u}}_I}$$

(51)

with

$${\mathbf{u}}_I={\left[w_I\hspace{0.17em}\hspace{0.17em}{\varphi }_I\right]}^T$$

(52)

where NP = p + 1 indicates the number of control points in each element, $R_I$ and ${\mathbf{u}}_I$ denote the NURBS basis function and the displacement vector at control point I, respectively, and $w_I$ and ${\varphi }_I$ are the deflection and angle of rotation of control point I, respectively.

Considering the strain expressions from Equations (12)–(14) and (15), the matrix form of the strain expressions can be described as

$${\epsilon }_{xx}=-z\frac{\partial \phi }{\partial x}={\mathit{C}}_1{\mathit{\epsilon }}_1$$

(53)

$${\epsilon }_{xz}=\frac{1}{2}\left(\frac{\partial w}{\partial x}-\phi \right)={\mathit{C}}_2{\mathit{\epsilon }}_2$$

(54)

$${\epsilon }_{yy}={\epsilon }_{zz}={\epsilon }_{xy}={\epsilon }_{yz}=0$$

(55)

with

$${\mathit{C}}_1=\left[1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}z\right]$$

(56)

$${\mathit{C}}_2=\left[1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1\right]$$

(57)

$${\mathit{\epsilon }}_1=\left[\begin{array}{c}0\\ -\frac{\partial \phi }{\partial x}\end{array}\right]$$

(58)

$${\mathit{\epsilon }}_2=\left[\begin{array}{c}\hspace{0.17em}\hspace{0.17em}\frac{1}{2}\frac{\partial w}{\partial x}\\ -\frac{1}{2}\phi \end{array}\right]$$

(59)

Using Equations (12)–(14), (16) and (17), the components of the rotation vector and the curvature tensor are obtained as

$$\mathit{\theta }=\left[\begin{array}{c}{\theta }_x\\ {\theta }_y\\ {\theta }_z\end{array}\right]=\left[\begin{array}{c}0\\ -\frac{1}{2}\left(\frac{\partial w}{\partial x}+\phi \right)\\ 0\end{array}\right]$$

(60)

$${\chi }_{xy}=-\frac{1}{4}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right)=C_2\mathit{\chi }$$

(61)

$${\chi }_{xx}={\chi }_{yy}={\chi }_{zz}={\chi }_{xz}={\chi }_{yz}=0$$

(62)

with

$$\mathit{\chi }=\left[\begin{array}{c}-\frac{1}{4}\frac{{\partial }^{2}w}{\partial x^2}\\ -\frac{1}{4}\frac{\partial \phi }{\partial x}\end{array}\right]$$

(63)

From Equations (12)–(14) and (18), the matrix form of the symmetric part of the second-order displacement gradient tensor can be written as

$${\eta }_{xxz}^s=\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}-\frac{2}{3}\frac{\partial \phi }{\partial x}={\mathit{C}}_2{\mathit{\eta }}_1$$

(64)

$${\mathit{\eta }}_1=\left[\begin{array}{c}\hspace{0.17em}\hspace{0.17em}\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}\\ -\frac{2}{3}\frac{\partial \phi }{\partial x}\end{array}\right]$$

(65)

For each element, by substituting Equations (51) and (52) into Equations (56)–(65) we obtain

$${\mathit{\epsilon }}_1={\displaystyle \sum _{I=1}^{NP}{\mathit{B}}_I^1}{\mathit{u}}_I$$

(66)

$${\mathit{\epsilon }}_2={\displaystyle \sum _{I=1}^{NP}{\mathit{B}}_I^2}{\mathit{u}}_I$$

(67)

$$\mathit{\chi }={\displaystyle \sum _{I=1}^{NP}{\mathit{B}}_I^m}{\mathit{u}}_I$$

(68)

$${\mathit{\eta }}_1={\displaystyle \sum _{I=1}^{NP}{\mathit{B}}_I^{s1}}{\mathit{u}}_I$$

(69)

at control point I

$${\mathit{B}}_I^1=\left[\begin{array}{cc}0& 0\\ 0& -\frac{d{R}_I}{dx}\end{array}\right]$$

(70)

$${\mathit{B}}_I^2=\left[\begin{array}{cc}\frac{1}{2}\frac{d{R}_I}{dx}& 0\\ 0& -\frac{1}{2}{R}_I\end{array}\right]$$

(71)

$${\mathit{B}}_I^m=\left[\begin{array}{cc}-\frac{1}{4}\frac{d^2{R}_I}{d{x}^2}& 0\\ 0& -\frac{1}{4}\frac{d{R}_I}{dx}\end{array}\right]$$

(72)

$${\mathit{B}}_I^{s1}=\left[\begin{array}{cc}\frac{1}{3}\frac{d^2{R}_I}{d{x}^2}& 0\\ 0& -\frac{2}{3}\frac{d{R}_I}{dx}\end{array}\right]$$

(73)

4.3. Static Bending

Considering the static bending problem, the body forces are ignored and all of the time derivatives are assumed to be zero. Using Equations (4), (12)–(21) and (29), the equation of motion can be rewritten as

$$\begin{array}{l}{\displaystyle {\int }_V\left[\delta {\epsilon }_{xx}\frac{E(1-v)}{(1+v)(1-2v)}{\epsilon }_{xx}+\delta {\epsilon }_{xz}\frac{E}{(1+v)}{\epsilon }_{xz}+\delta {\chi }_{xy}\frac{E{l}_m^2}{(1+v)}{\chi }_{xy}+\delta {\eta }_{xxz}^s\frac{E{l}_s^2}{(1+v)}{\eta }_{xxz}^s\right]}dV\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle {\int }_0^l\delta w{q}_{z}dx+\delta wP}\hspace{0.17em}\end{array}$$

(74)

where q_z represents the distributed load and P denotes the concentrated load. By inserting Equations (53)–(65) into Equation (74), we obtain

$$\begin{array}{l}{{\displaystyle {\int }_0^L\delta \left[\begin{array}{c}{\mathit{\epsilon }}_1\\ {\mathit{\epsilon }}_2\end{array}\right]}}^T\left[\begin{array}{cc}{\mathit{D}}_1& 0\\ 0& {\mathit{D}}_2\end{array}\right]\left[\begin{array}{c}{\mathit{\epsilon }}_1\\ {\mathit{\epsilon }}_2\end{array}\right]dx+{\displaystyle {\int }_0^L\delta {\mathit{\eta }}_1{}^T}{\mathit{D}}_{s1}{\mathit{\eta }}_{1}dx+{\displaystyle {\int }_0^L\delta {\mathit{\chi }}^T{\mathit{D}}_m\mathit{\chi }dx}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle {\int }_0^L\delta q_{z}wdx+\delta Pw}\end{array}$$

(75)

with

$${\mathit{D}}_1=\underset{A}{{\displaystyle \int }}(2\mu +\lambda )\left[\begin{array}{cc}0& 0\\ 0& z^2\end{array}\right]dA$$

(76)

$${\mathit{D}}_2=\underset{A}{{\displaystyle \int }}4\mu K_s\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]dA$$

(77)

$${\mathit{D}}_m=\underset{A}{{\displaystyle \int }}4\mu l_m^2\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]dA$$

(78)

$${\mathit{D}}_{s1}=\underset{A}{{\displaystyle \int }}6\mu l_s^2\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]dA$$

(79)

Considering the rectangular cross-section of a microbeam with width b and height h, Equations (76)–(79) can be rewritten as follows:

$${\mathit{D}}_1=\frac{E(1-v)}{(1+v)(1-2v)}\left[\begin{array}{cc}0& 0\\ 0& \frac{b{h}^3}{12}\end{array}\right]$$

(80)

$${\mathit{D}}_2=\frac{2Ebh{K}_s}{1+v}\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]$$

(81)

$${\mathit{D}}_m=\frac{2Ebh{l}_m^2}{1+v}\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]$$

(82)

$${\mathit{D}}_{s1}=\frac{3Ebh{l}_s^2}{1+v}\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]$$

(83)

From Equation (75), the equation of motion of static problem can be obtained as

$$\mathbf{K}\mathbf{d}=\mathbf{F}$$

(84)

where d is the displacement vector of the control points and the global stiffness matrix, K, and the load vector, F, can be obtained by assembling each element written as

$$\begin{array}{l}\mathbf{K}={\displaystyle \sum _{e=1}^{Ne}{\displaystyle \sum _{I=1}^{NP}{{\displaystyle {\int }_0^L\left[\begin{array}{c}{\mathit{B}}_I^1\\ {\mathit{B}}_I^2\end{array}\right]}}^T\left[\begin{array}{cc}{\mathit{D}}_1& 0\\ 0& {\mathit{D}}_2\end{array}\right]\left[\begin{array}{c}{\mathit{B}}_I^1\\ {\mathit{B}}_I^2\end{array}\right]dx}}+{\displaystyle \sum _{e=1}^{Ne}{\displaystyle \sum _{I=1}^{NP}{\displaystyle {\int }_0^L{({\mathit{B}}_I^{s1})}^T}{\mathit{D}}_{s1}{\mathit{B}}_I^{s1}dx}}\\ +{\displaystyle \sum _{e=1}^{Ne}{\displaystyle \sum _{I=1}^{NP}{\displaystyle {\int }_0^L{({\mathit{B}}_I^m)}^T{\mathit{D}}_m{\mathit{B}}_I^m}dx}}\end{array}$$

(85)

$$\mathbf{F}={\displaystyle \sum _{e=1}^{Ne}{\displaystyle \sum _{I=1}^{NP}{\displaystyle {\int }_0^L{\mathit{B}}_I^f{q}_{z}dx}}}+{\displaystyle \sum _{e=1}^{Ne}{\displaystyle \sum _{I=1}^{NP}{\mathit{B}}_I^{f}P}}$$

(86)

where Ne is the total number of elements and where

$${\mathit{B}}_I^f={\left[R_I\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\right]}^T$$

(87)

4.4. Free Vibration

For the free vibration analysis, all external forces are set to zero. Based on the variational method, from Equations (25), (30) and (74) the governing equation can be described as

$$\begin{array}{l}{\displaystyle {\int }_V\left[\delta {\epsilon }_{xx}\frac{E(1-v)}{(1+v)(1-2v)}{\epsilon }_{xx}+\delta {\epsilon }_{xz}\frac{E}{(1+v)}{\epsilon }_{xz}+\delta {\chi }_{xy}\frac{E{l}_m^2}{(1+v)}{\chi }_{xy}+\delta {\eta }_{xxz}^s\frac{E{l}_s^2}{(1+v)}{\eta }_{xxz}^s\right]}dV\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle {\int }_V\left[\delta \phi z^2\rho \ddot{\phi }+\delta w\rho \ddot{w}+\delta \left(\frac{\partial \phi }{\partial x}\right)z^2\rho l_v^2\frac{\partial \ddot{\phi }}{\partial x}+\delta \phi \rho l_v^2\ddot{\phi }+\delta \left(\frac{\partial w}{\partial x}\right)\rho l_v^2\frac{\partial \ddot{w}}{\partial x}\right]}\hspace{0.17em}dV\end{array}$$

(88)

Based on D’Alembert’s principle and substituting Equations (51)–(73) into (88), we obtain

$$\left(\mathbf{K}-{\omega }^2\mathit{M}\right)\mathit{d}=0$$

(89)

in which ω denotes the natural frequency and the global mass matrix M, obtained by assembling each element, can be described as

$$\mathit{M}={\displaystyle \sum _{e=1}^{Ne}{\displaystyle \sum _{I=1}^{NP}{\displaystyle {\int }_0^L\left({\left[\begin{array}{cc}{R}_I& 0\\ 0& R_I\end{array}\right]}^T\left[\begin{array}{cc}\rho A& 0\\ 0& I\rho \end{array}\right]\left[\begin{array}{cc}{R}_I& 0\\ 0& R_I\end{array}\right]+{(R_I)}^T\left[\begin{array}{ccc}\rho A{l}_{{}^v}^2& 0& 0\\ 0& \rho I{l}_v^2& 0\\ 0& 0& \rho A{l}_v^2\end{array}\right]{\mathit{R}}_I\right)dx}}}$$

(90)

at control point I

$${\mathit{R}}_I=\left[\begin{array}{cc}{R}_{I,x}& 0\\ 0& R_{I,x}\\ 0& R_I\end{array}\right]$$

(91)

5. Examples and Discussions

In order to certify the present non-classical Timoshenko–Ehrenfest beam model, several examples of static bending and free vibration problems (see Figure 2a) were solved by both analytical solution and the isogeometric analysis approach.

Using Equations (34)–(36), the simply-supported boundary conditions can be described as follows:

$${w|}_{x=0}={w|}_{x=L}=0$$

(92)

$$M_x+\frac{1}{2}{Y}_x+2Q_x-l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}=0, \mathrm{at} x=0 \mathrm{and} x=L$$

(93)

$$\frac{1}{2}{Y}_x-Q_x=0, \mathrm{at} x=0 \mathrm{and} x=L$$

(94)

$${\frac{{\partial }^{2}w}{\partial x^2}|}_{x=0}={\frac{{\partial }^{2}w}{\partial x^2}|}_{x=L}=0$$

(95)

which is applied through isogeometric analysis unless otherwise specified.

Using Equations (37)–(40), Equations (92) and (94) can be rewritten as

$$\begin{array}{l}4l_s^2\mu A\left(\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}-\frac{2}{3}\frac{\partial \phi }{\partial x}\right)-\frac{E(1-v)I}{(1+v)(1-2v)}\frac{\partial \phi }{\partial x}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{1}{4}\mu l_m^{2}A\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right)-l_v^2{m}_2\frac{{\partial }^3\phi }{\partial x\partial t^2}\hspace{0.17em}=0\end{array}$$

(96)

$$\frac{1}{4}\mu l_m^{2}A\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right)+2l_s^2\mu A\left(\frac{1}{3}\frac{{\partial }^{2}w}{\partial x^2}-\frac{2}{3}\frac{\partial \phi }{\partial x}\right)=0$$

(97)

5.1. Static Bending Problem

Considering the deformation as foreign to the time t and neglecting the time-dependent partial derivatives, the static bending boundary value problem of simply-supported microbeam (as shown in Figure 2a) is described by Equations (42) and (43) and (92)–(97), with φ = φ(x) and w = w(x). In addition, the body couple c is set to zero.

The deformations w(x) and φ(x) are described by the following Fourier series solutions:

$$w(x)={\displaystyle \sum _{n=1}^{\infty }{W}_n}\sin \left(\frac{n\pi x}{L}\right)\hspace{0.17em}$$

(98)

$$\phi (x)={\displaystyle \sum _{n=1}^{\infty }{\Phi }_n}\cos \left(\frac{n\pi x}{L}\right)$$

(99)

where $W_n$ and ${\Phi }_n$ are Fourier coefficients defined for each n. Obviously, the formulae w(x) and φ(x) from Equation (62) satisfy the boundary conditions in Equations (92)–(97) for any $W_n$ and any ${\Phi }_n$.

According to Equations (42) and (92)–(95), the external load, q(x), is described by Fourier series as

$$q(x)={\displaystyle \sum _{n=1}^{\infty }{Q}_n}\sin \left(\frac{n\pi x}{L}\right)$$

(100)

where the Fourier coefficient, $Q_n$, can be described as

$$Q_n=\frac{2}{L}{\displaystyle {\int }_0^{L}q(x)\sin \left(\frac{n\pi x}{L}\right)}dx$$

(101)

In the static bending problem, we consider q(x) = Pδ(x − L/2) (see Figure 2a), with δ(.) denoting the Dirac function. By substituting q(x) in this paragraph into Equation (101), we obtain

$$Q_n=\frac{2}{L}P\sin \left(\frac{n\pi }{2}\right)$$

(102)

Based on the above analysis, by inserting Equations (98)–(100) into Equations (42) and (43) we obtain

$$\begin{array}{l}\left[\frac{2}{3}\mu l_s^{2}A{\left(\frac{n\pi }{L}\right)}^4+\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^4+K_s\mu A{\left(\frac{n\pi }{L}\right)}^2\right]W_n+[\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^3-K_s\mu A\left(\frac{n\pi }{L}\right)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{4}{3}\mu l_s^{2}A{\left(\frac{n\pi }{L}\right)}^3]{\Phi }_n=Q_n\end{array}$$

(103)

$$\begin{array}{l}\left[\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^3-K_s\mu A\left(\frac{n\pi }{L}\right)-\frac{4}{3}\mu l_s^{2}A{\left(\frac{n\pi }{L}\right)}^3\right]W_n+[\frac{8}{3}\mu l_s^{2}A{\left(\frac{n\pi }{L}\right)}^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{E(1-v)I}{(1+v)(1-2v)}{\left(\frac{n\pi }{L}\right)}^2+K_s\mu A+\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^2]{\Phi }_n=0\end{array}$$

(104)

By solving the two existing non-vanishing governing equations, we can acquire

$${\Phi }_n=-\frac{Q_{{}_n}\left(-\frac{4}{3}{l}_s^2-K_s{\left(\frac{L}{n\pi }\right)}^2+\frac{1}{4}{l}_m^2\right)L^3}{{\left(n\pi \right)}^3\left(\frac{3}{2}{l}_s^2{l}_m^2\mu A+C_1+C_2\right)}$$

(105)

$$W_n=\frac{Q_n\left(\frac{8}{3}\mu l_s^{2}A+\frac{E(1-v)I}{(1+v)(1-v)}+K_s\mu A{\left(\frac{L}{n\pi }\right)}^2+\frac{1}{4}\mu l_m^{2}A\right)L^4}{{\left(n\pi \right)}^4\mu A\left(\frac{3}{2}{l}_s^2{l}_m^2\mu A+C_1+C_2\right)}$$

(106)

with

$$C_1=\left(\frac{2}{3}{l}_s^2+l_m^2\right)\mu A{K}_s{\left(\frac{L}{n\pi }\right)}^2$$

(107)

$$C_2=\left[\frac{2}{3}{l}_s^2+\frac{1}{4}{l}_m^2+K_s{\left(\frac{L}{n\pi }\right)}^2\right]\frac{EI(1-v)}{(1+v)(1-2v)}$$

(108)

Then, the analytical expressions of w(x) and φ(x) can be acquired by substituting Equations (105) and (106) into Equations (98) and (99), respectively. When the material parameter l_s is set to be zero, Equations (105)–(108) degenerates to

$${\Phi }_n=-\frac{Q_n\left(n^2{\pi }^2{l}_m^2-4K_s{L}^2\right)L^3}{{\left(n\pi \right)}^3\left(4l_m^2\mu A{K}_s{L}^2+\frac{EI(1-v)}{(1+v)(1-2v)}\left(l_m^2{n}^2{\pi }^2+4K_s{L}^2\right)\right)}$$

(109)

$$W_n=\frac{Q_n\left(\frac{4n^2{\pi }^{2}EI(1-v)}{(1+v)(1-2v)}+4K_s\mu A{L}^2+l_m^{2}A{n}^2{\pi }^2\right)L^4}{{\left(n\pi \right)}^4\mu A\left(4l_m^2\mu A{K}_s{L}^2+\frac{EI\left(1-v\right)}{\left(1+v\right)\left(1-2v\right)}\left(l_m^2{n}^2{\pi }^2+4K_s{L}^2\right)\right)}$$

(110)

which are identical to those rotation and deflection expressions provided by the modified couple stress theory [38]. Both the analytic solutions and isogeometric analysis solutions of the change of deflection are presented in Figure 3 for comparison. In the isogeometric analysis, the NURBS basis functions with degree of p = 3 and 30 control points are applied, which can be obtained from the following case of Figure 4 and Figure 5. In the current case, the concentrated load P = 100 µN is applied to a simply-supported microbeam (see Figure 2a). The configuration of the microbeam is defined by letting L = 20 h, b = 2 h, and the thickness h = 17.6 µm. The epoxy is considered in this work, and the parameters of the beam material involved in the numerical study examples study are E = 1.44 GPa, v = 0.38. The material scale parameter of the couple stress effect is l_m = 17.6 µm [1,36]. As shown in Figure 3, all the deflections obtained by the presented non-classical Timoshenko-Ehrenfest beam model using IGA with different l_s (l_s = 0, 0.3l_m, 0.6l_m, 1.2l_m) and by its classical beam model are in remarkable agreement with the analytic solution obtained from Equations (105)–(108). It can be seen from Figure 3 that the deflection, w/h, obtained by the present model based on RSGT is always lower than that predicted by the classical model, which shows the behavior of the size effect. In addition, the deflection acquired by the proposed model decreases as the material parameter l_s increases.

To verify the convergence of the IGA based on RSGT, the relative error in L²-norm with respect to element size of a simply-supported beam for different orders p = 2, 3, 4, 5, 6 with thickness values h = 20 μm, 40 μm are studied here, as shown in Figure 4; the analytical solution obtained from Equations (98) and (99) is incorporated into this convergence analysis for comparison. From Figure 4, it can be observed that the acquired convergence results improve as the element size decreases, and the higher orders (p = 3, 4, 5, 6) of IGA have a better convergence rate than p = 2, which is in line with the expectation that a higher order results in better precision and convergence rate.

Figure 5 shows the maximum deflection w/h obtained from the IGA with different degrees p = 2, 3, 4 and number of control points for thickness h = 20 μm, 40 μm. It can be seen from Figure 5 that the numerical results for the maximum deflection w/h predicted by IGA for order p = 3 and p = 4 with the number of control point 25 are essentially in agreement with the analytical solution. According to the above conclusion, the cubic NURBS basis function with control point 30 essentially fulfills the convergence requirements, and is adopted in the rest of this work. In addition, the four Gaussian integration points are adopted in the NURBS element except for other specified special circumstances. In addition, the assumption in the expansion of Equations (98) and (99) that the numerical value n is 30 is sufficient for the convergent analysis, as shown in Figure 3.

To verify the accuracy of the proposed model, a comparative study of the variations of deflection between Euler–Bernoulli theory (EBT) as per Zhang [44] and Timoshenko–Ehrenfest beam theory (TBT) as per our proposed model were considered, as shown in Figure 6. The configuration and material properties of the microbeam are the same as mentioned above, with h = 17.6 μm. Figure 6 shows that the variations of the deflection obtained by Timoshenko–Ehrenfest beam model using IGA are in remarkable agreement with the analytic solution obtained from Equations (105)–(108), and the deflection obtained by EBT is always lower than that obtained by TBT. In addition, it can be observed from Figure 6 that the deviation of the results obtained by EBT and TBT is negligible for a slender beam (here, h = 17.6 µm and L = 20 h).

Figure 7 presents the variations of deflection w/h obtained by the reformulated strain gradient theory (RSGT) with l_s = 1.2 l_m, modified couple stress theory (MCST) and classical Timoshenko–Ehrenfest beam theory (Classical) with h = l, 2l (l_m = l = 17.6 µm). From Figure 7, it can be clearly observed that not only is the deflection obtained by RSGT always lower than that obtained by classical theory, it is lower than that obtained by MCST for all cases, which indicates that in the non-classical Timoshenko–Ehrenfest microbeam model based on RSGT, the bending rigidity appears to increase as the material scale parameter l_s is introduced.

Figure 8 shows influence of different thickness values (with h = l, 2l, 4l and 8l) on deflection as predicted by the current Timoshenko–Ehrenfest beam model with three different theories (including RSGT with l_s = 1.2l_m, MCST, and Classical theory). It can be clearly observed from Figure 8 that the deflection predicted by the current Timoshenko–Ehrenfest beam model based on RSGT and MCST is invariably lower than that predicted by classical Timoshenko–Ehrenfest beam theory in any instance. Figure 8 further shows that the distinction between the deflection obtained with the non-classical Timoshenko–Ehrenfest beam theory and with classical Timoshenko-Ehrenfest beam theory is significant only with a very small beam thickness (here, with h < 2l). However, the distinction narrows as the thickness value becomes greater (here, up to 4l). This observation indicates that the size effect is remarkable as thickness decreases, which is in agreement with the general tendency observed experimentally [2].

To certify the applicability of the proposed isogeometric analysis approach to a common problem, a static bending analysis of a cantilever beam subject to concentrated load P = 0.01N and a clamped–clamped beam under a concentrated load P = 0.01N (see Figure 2) was investigated. The geometry of the microbeams was defined as width b = 2 h, length L = 20 h and thickness h = 17.6 µm, while beam materials involved in this particular numerical example were the same as those adopted in Section 5.1. Figure 9 and Figure 10 show the results of the deflection w/h predicted by IGA based on RSGT (with l_s = 0.3l_m, 0.6l_m, 1.2l_m), MCST, and classical Timoshenko–Ehrenfest beam theory (Classical) for a cantilever beam and clamped–clamped beam with a beam thickness value h (h = l, 2l, 4l and l_m = l = 17.6 µm). It can be seen from Figure 9 and Figure 10 that the results obtained with the IGA based on RSGT and MCST is always smaller than those obtained by classical theory for all cases, and the deflection predicted by the current numerical approach decreases with the increase of the strain gradient, $l_s$. It can be further concluded from Figure 9 and Figure 10 that the differences between the deflection values predicted by current non-classical Timoshenko–Ehrenfest beam theory and those obtained by the classical Timoshenko–Ehrenfest beam theory are significant, as the beam thickness, h, tends to be small.

5.2. Free Vibration Problem

This section investigates the free vibration problem of a simply-supported microbeam. The boundary conditions are provided in Equations (92) and (95), with the external forces disappeared (i.e., q = 0 and c = 0; $\overline{N}=0$; and $\overline{M}=0$ and $\overline{H}=0$). Conforming to the static bending analysis mentioned above, the following Fourier series expansions for w(x,t) and φ(x,t) are applied:

$$w(x,t)={\displaystyle \sum _{n=1}^{\infty }{W}_n^V}\sin \left(\frac{n\pi x}{L}\right)e^{i{\omega }_{n}t}$$

(111)

$$\phi (x,t)={\displaystyle \sum _{n=1}^{\infty }{\Phi }_n^V}\cos \left(\frac{n\pi x}{L}\right)e^{i{\omega }_{n}t}$$

(112)

where $W_n^V$ and ${\Phi }_n^V$ are Fourier coefficients, ω_n denotes the vibration frequency, and i is a general imaginary number observing i² = −1. Obviously, the Fourier series expansions in Equation (69) observe the boundary conditions in Equations (92) and (95) for any $W_n^V$ and ${\Phi }_n^V$. Substituting Equation (111) into Equations (42) and (43) provides

$$\begin{array}{l}\left[m_0\left(-{\omega }_n^2-l_v^2{\omega }_n^2{\left(\frac{n\pi }{L}\right)}^2\right)+K_s\mu A{\left(\frac{n\pi }{L}\right)}^2+\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^4+\frac{2}{3}{l}_s^2\mu A{\left(\frac{n\pi }{L}\right)}^4\right]W_n^V\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left[-K_s\mu A\left(\frac{n\pi }{L}\right)+\frac{1}{4}\mu A{l}_m^2{\left(\frac{n\pi }{L}\right)}^3-\frac{4}{3}{l}_s^2\mu A{\left(\frac{n\pi }{L}\right)}^3\right]{\Phi }_n^V=0\end{array}$$

(113)

$$\begin{array}{l}\left[-K_s\mu A\left(\frac{n\pi }{L}\right)+\frac{1}{4}\mu A{l}_m^2{\left(\frac{n\pi }{L}\right)}^3-\frac{4}{3}{l}_s^2\mu A{\left(\frac{n\pi }{L}\right)}^3\right]W_n^V+[-m_2\left({\omega }_n^2+l_v^2{\omega }_n^2{\left(\frac{n\pi }{L}\right)}^2\right)-m_0{l}_v^2{\omega }_n^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{EI(1-v)}{(1+v)(1-2v)}{\left(\frac{n\pi }{L}\right)}^2+K_s\mu A+\frac{1}{4}\mu A{l}_m^2{\left(\frac{n\pi }{L}\right)}^2+\frac{8}{3}\mu A{l}_s^2{\left(\frac{n\pi }{L}\right)}^2]{\Phi }_n^V=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}$$

(114)

Using Equations (113) and (114) with a non-trivial solution of $W_n^V\ne 0$ and ${\Phi }_n^V\ne 0$, we can obtain

$$S_0{\omega }_n^4-(S_1+S_2){\omega }_n^2+S_3=0$$

(115)

with

$$S_0=\left[m_0+m_0{l}_v^2{\left(\frac{n\pi }{L}\right)}^2\right]\left[m_2+m_2{l}_v^2{\left(\frac{n\pi }{L}\right)}^2+m_0{l}_v^2\right]$$

(116)

$$S_1=\left[m_0+m_0{l}_v^2{\left(\frac{n\pi }{L}\right)}^2\right][\frac{EI(1-v)}{(1+v)(1-2v)}{\left(\frac{n\pi }{L}\right)}^2+K_s\mu A+\frac{1}{4}\mu A{l}_m^2{\left(\frac{n\pi }{L}\right)}^2+\frac{8}{3}\mu A{l}_s^2{\left(\frac{n\pi }{L}\right)}^2]$$

(117)

$$S_2=\left[m_2+m_2{l}_v^2{\left(\frac{n\pi }{L}\right)}^2+m_0{l}_v^2\right]\left[K_s\mu A{\left(\frac{n\pi }{L}\right)}^2+\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^4+\frac{2}{3}{l}_s^2\mu A{\left(\frac{n\pi }{L}\right)}^4\right]$$

(118)

$$\begin{array}{l}{S}_3=\mu A{\left(\frac{n\pi }{L}\right)}^6[\left(\frac{K_s{L}^2}{n^2{\pi }^2}+\frac{1}{4}{l}_m^2+\frac{2}{3}{l}_s^2\right)\frac{EI(1-v)}{(1+v)(1-2v)}+l_m^2{K}_s\mu A{\left(\frac{L}{n\pi }\right)}^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{2}{3}{l}_s^2{K}_s\mu A{\left(\frac{L}{n\pi }\right)}^2\hspace{0.17em}+\frac{3}{2}{l}_s^2{l}_m^2\mu A]\hspace{0.17em}\end{array}$$

(119)

The solution of Equation (115) with respect to the quadratic equation of ${\omega }_n^2$ can be obtained as

$${\omega }_n^2=\frac{S_1+S_2-\sqrt{{(S_1+S_2)}^2-4S_0{S}_3}}{2S_0}$$

(120)

When the material scale parameter l_s and velocity gradient coefficient l_v are assumed to be zero, Equations (116)–(119) can be further degenerated to

$$S_0=m_0{m}_2$$

(121)

$$S_1=m_0[\frac{EI(1-v)}{(1+v)(1-2v)}{\left(\frac{n\pi }{L}\right)}^2+K_s\mu A+\frac{1}{4}\mu A{l}_m^2{\left(\frac{n\pi }{L}\right)}^2]$$

(122)

$$S_2=m_2\left[K_s\mu A{\left(\frac{n\pi }{L}\right)}^2+\frac{1}{4}\mu l_m^{2}A{\left(\frac{n\pi }{L}\right)}^4\right]$$

(123)

$$S_3={\left(\frac{n\pi }{L}\right)}^4[\frac{E{K}_s\mu (1-v)IA}{(1+v)(1-2v)}+l_m^2{K}_s{\mu }^2{A}^2+\frac{1}{4}\frac{E{l}_m^2\mu (1-v)IA}{(1+v)(1-2v)}{\left(\frac{n\pi }{L}\right)}^2]$$

(124)

Substituting Equation (121) into Equation (120) then obtains the expressions of the first natural frequency provided by modified couple stress theory [38]. To certify the accuracy of the first natural frequency provided by IGA, the results of the first natural frequency obtained using analytic solutions, from Equation (120) with n = 1, and the numerical solutions obtained from IGA are presented in Figure 11 for comparison. Figure 11 shows the variations of the first natural frequency predicted by the non-classical Timoshenko–Ehrenfest beam analytical solutions based on RSGT (with l_s = 1.2l_m) for different velocity gradient coefficients l_v [22,25,44], modified couple stress theory (with l_v = l_s = 0), classical Timoshenko–Ehrenfest beam theory (with l_v = l_s = l_m = 0), and the numerical result provided by IGA. The material parameter and configuration (with L = 20 h, b = 2 h, and h = 17.6 µm) of the simply-supported beam (epoxy) considered here is the same as that applied in Section 5.1. The material density is ρ = 1.22 × 10³ Kg/m³.

In Figure 11, it can be clearly seen that the analytical result obtained by the current Timoshenko–Ehrenfest beam model agrees well with the numerical result obtained by IGA, and the natural frequency diminishes with the increase of the velocity gradient coefficient, l_v. The distinction between the natural frequency obtained with the current Timoshenko–Ehrenfest beam model based on non-classical theories (RSGT and MCST with l_v) and classical Timoshenko–Ehrenfest beam theory are significant when the thickness value is small (h < 2l_m); the differences narrow as h becomes greater. This certifies the micro-dependent size effect for free vibration. In addition, the natural frequency predicted by both the Timoshenko–Ehrenfest beam model based on RSGT (here, with l_v < 12l_m) and MCST are always higher than that obtained using classical Timoshenko–Ehrenfest beam theory. However, with a velocity gradient coefficient up to 12l_m (here, l_v ≥ 12l_m), the natural frequency obtained using the Timoshenko–Ehrenfest beam model based on RSGT is lower than that obtained using classical Timoshenko–Ehrenfest beam theory. This roughly reflects the effect of three additional scale parameters (l_s, l_m, and l_v). To further illustrate the effect of material scale parameters, the results shown in Figure 12 are applied in following numerical example.

Figure 12 shows the effect of three material scale parameters (l_s, l_m, and l_v, denoting the strain gradient, couple stress, and velocity gradient effect, respectively) on the natural frequency obtained by the current Timoshenko–Ehrenfest beam model. From Figure 12, it can be observed that the natural frequency increases as l_s becomes greater, while the natural frequency diminishes as l_v increases. This indicates that the velocity gradient coefficient decreases with the natural frequency; however, the strain gradient, l_s, and couple stress, l_m, cause the natural frequency to increase; these influences are only distinct when the thickness value, h, is small (here, h ≤ 2l_m).

6. Conclusions

This work proposes a new non-classical Timoshenko–Ehrenfest beam model based on a reformulated strain gradient elasticity theory. The couple stress, strain gradient, and velocity gradient effects are considered in the reformulated strain gradient elasticity theory by only one material scale parameter for each. The static bending and free vibration problems are solved for Timoshenko–Ehrenfest beams using both analytical solutions and the isogeometric analysis approach with high-order continuity non-uniform rational B-spline basis functions, which effectively fulfills the higher derivatives requirement in strain gradient theory. Modified couple stress theory and classical elasticity theory can be considered as special cases in the proposed non-classical Timoshenko–Ehrenfest beam formulation by suppressing the corresponding strain gradient, couple stress and velocity gradient effects, respectively. In order to demonstrate the proposed non-classical Timoshenko–Ehrenfest beam model and verify the accuracy and efficiency of the proposed isogeometric approach, several examples of static bending and free vibration problems were studied and compared with each other.

It can be concluded from the numerical examples that the strain gradient and couple stress effect in the proposed non-classical model result in decreased deflection and an increased natural frequency. On the other hand, the microstructure-dependent size effects on the natural frequency are significant when the beam thickness is very small, and the presence of the velocity gradient effect leads to a decrease in the natural frequency. The numerical results show that the observed size effects in the proposed non-classical Timoshenko–Ehrenfest beam model match very well with observations of the same phenomenon in experiments and analytical solutions.