Nonlinear free vibration of size-dependent microbeams with nonlinear elasticity under various boundary conditions

In this study, nonlinear couple stress–strain constitutive relationships in the modified couple stress theory (MCST) are derived on the basis of previous classical stress–strain constitutive relationships of nonlinear elasticity materials. Hamilton’s principle is employed to obtain higher-ordernonlineargoverningequationswithintheframeworkoftheupdatedMCST,vonKármángeometricnonlinearityandBernoulli–Eulerbeamtheory.Thesemathematicalformulationsaresolvednumericallybythedifferentialquadraturemethodtogetherwithaniterativealgorithmtodeterminethenonlineardynamicfeaturesofmicrobeamswithfourgroupsofboundaryconditions.Adetailedparametricstudyisconductedtoanalyzetheinfluencesofnonlinearelasticitypropertiesonthenonlinearfreevibrationcharacteristicsofthemicrobeams.Resultsshowthatthesemicrobeamsexhibitingnonlinearcoupleconstitutiverelationshipshavelowerfrequenciesthantheirapproximatelysimplifiedlinearcoupleconstitutiverelationships.Inaddition,thefrequenciesofmicrobeamswithnonlinearelasticitypropertiesdecreaseasthevibrationamplitudeincreases.

In recent years, studies have applied the nano-and microbeams on linear and nonlinear free vibrations to predict dynamic features such as natural frequencies and mode shapes. Kong et al. [9] established a model to solve analytically the dynamic problems of Bernoulli-Euler beams on the basis of modified couple stress theory (MCST). They concluded that linear natural frequencies are size dependent and that the difference of natural frequencies between the new and the classical beam model is very significant. On the basis of the sinusoidal beam theory and the modified strain gradient theory, Wang et al. [10] studied size-dependent bending and vibration of microbeams made of porous metal foams. Wang and Wang [11] first presented an exact linear free vibration solution for exponentially tapered cantilever with tip mass and obtained the exact results, i.e. linear frequencies and mode shapes. Thereafter, Wang [12] investigated the linear free vibration of a tapered cantilever of constant thickness and linearly tapered width. Nourbakhsh et al. [13] investigated the problems of nonlinear free and forced vibrations of microscale simply supported beams with initial lateral displacement. For free vibrations, they concluded that the nonlinear frequencies are much higher than the linear ones. On the basis of Eringen's nonlocal elasticity theory and von Kármán's geometric nonlinearity, Ke et al. [14] studied the nonlinear free vibration of embedded double-walled carbon nanotubes (DWNTs) and investigated the influences of nonlocal parameter, tube length, spring constant and end supports on the nonlinear free vibration characteristics of DWNTs. They found that both linear frequency and nonlinear frequency ratios decrease as the length of DWNTs increases. Subsequently, Ke et al. [15] investigated the nonlinear free vibration of microbeams made from functionally graded materials (FGMs) on the basis of the MCST and von Kármán geometric nonlinearity. Their findings showed that linear and nonlinear frequencies increase significantly when the thickness of the FGM microbeam is comparable to the material length scale parameter. Asghari et al. [16] introduced a size-dependent nonlinear Timoshenko microbeam model on the basis of strain gradient theory to investigate nonlinear static and free vibration behaviors of hinged-hinged beams. Some methods are applicable for the nonlinear analysis of structures, such as the multiple-scale method [17,18], the harmonic balance method [19], and the differential quadrature method (DQM) [14,15,20,21].

THE MCST IN NONLINE AR EL A STICIT Y 2.1 The MCST in linear elasticity
On the basis of the MCST given by Yang et al. in 2002 [24], the strain energy density φ(ε ij , η ij ) for an arbitrary volume element is given by where σ ij , ε ij , m ij and η ij represent the stress tensor, strain tensor, deviatoric part of the couple stress tensor and symmetric curvature tensor, respectively. They are defined as follows: where γ i is the component of the rotation vector; l is a material length scale parameter; μ is Poisson's ratio; and λ, E and G are Lamé constants.

Shear stress-strain constitutive relationships in nonlinear elasticity
As stated before, Peng et al. [36] utilized the following third-order power function to fit the experimentally determined nonlinear normal stress-strain constitutive relationship curve of isotropic nonlinear elasticity material Al-1% Si [33]: where E 1 and E 2 are two generalized Young's modulus coefficients. Moreover, Eq. (2a) could be rewritten as where and σ α xx (α = L, N) represent linear and nonlinear components, respectively. Taking the first derivative of Eq. (2a) with respect to ε xx , we could get its differential equation written as where E xx is the generalized nonlinear elastic modulus, i.e. the gradient of Eq. (2a) at arbitrary strain ε xx , along the x direction and its expression is given as Derived from the nonlinear classical stress-strain constitutive relationship (2a), differential equation (3a) exhibits linear relationship between the differentials of the classical normal stress dσ xx and strain dε xx at arbitrary strain ε xx , so it follows Hooke's law. For the isotropic nonlinear elasticity materials, we could also have the following similar differential equations like Eqs (3a) and (3b) along the other two directions by replacing their double subscript xx with ii (i = y, z): where E ii is the generalized nonlinear elastic modulus along the i direction and is given as Now consider a cube element abcda b c d of isotropic nonlinear elasticity materials and it is only subjected to differential of plane stresses dσ βϕ (β, ϕ = x, y) as shown in Fig. 1a, where the three axes of Cartesian coordinate xyz are vertical to three planes of the cube element. Referring to the effect of the Poisson's ratio, the differential of strain dε xx in the x direction due to the differential of stress dσ xx is equal to dσ xx /E xx and in the y direction due to the differential of stress dσ yy is −μ dσ yy /E yy . Thus, the differential of the resultant strain in the x direction is In a similar manner, the differential of the strain dε yy in the y direction is obtained as However, for the element of isotropic nonlinear elasticity materials subjected to the differential of plane shear stresses dσ xy , we assume that the relationship between the differentials of the shear stress dσ xy and strain dε xy follows Hooke's law and the differential of shear strain dε xy is given as where G xy is the generalized shear modulus of the shear stress-strain constitutive relationships for the nonlinear elasticity materials and is unknown. Equations (5a)-(5c) could be rewritten for the differential of the stresses in terms of the differential of the strains as As we know, general equations of plane-stress transformation and plane-strain transformation are obtained on the basis of linear elasticity. However, Eqs (6a)-(6c) are the linear relationships between differentials of the stresses dσ αβ and strains dε αβ (α, β = x, y), and they follow Hooke's law. We cut a triangular subelement ebf along the z direction from the element abcda b c d and it is shown in Fig. 1b. Referring to general equations of plane-stress transformation in Mechanics of Materials [39], we could obtain general equations of differential-of-plane-stress transformation as follows: where θ is the angle between the normal directionx of the plane ef and the x direction. Similarly, referring to general equations of plane-strain transformation, general equations of differential-of-plane-strain transformation of the triangular subelement ebf could be obtained directly as On the basis of general equations of differential-of-plane-stress transformation, Eqs (7a)-(7c), the cube element abcd in Fig. 1c, which is only subjected to differentials of pure shear stresses dσ xy and dσ yx , is equal to the cube element a 1 b 1 c 1 d 1 subjected to differentials of principal stresses dσ 45 • xx and dσ 45 • yŷ at θ = 45°. Meanwhile, by setting θ = 45°in Eqs (7a)-(7c) and (8a)-(8c), one has the following differential equations of differential-of-plane-stress transformation and differential-of-plane-strain transformation: Integrating Eq. (10a) with respect to ε xy from zero to ε xy under the condition that ε 45 • xx = 0 at ε xy = 0, one obtains On the basis of Eq. (6a), differential of principal stress dσ 45 • xx is obtained as By substituting Eqs (9a), (10a), (10b) and (10d) into Eq. (11a), one obtains Integrating the above equation with respect to ε xy from zero to ε xy under the condition that σ xy = 0 at ε xy = 0, the nonlinear classical shear stress-strain constitutive relationship is obtained as For the isotropic nonlinear elasticity materials, the nonlinear classical shear stress-strain constitutive relationships could be obtained by repeating the above process and their expressions are given as

Couple stress-strain constitutive relationships in nonlinear elasticity
We have obtained the nonlinear classical shear stress-strain constitutive relationships Eq. (12) for the nonlinear elasticity materials. Next, we are required to determine the expressions of the couple stress-strain constitutive relationships for the nonlinear elasticity materials.
First, let us obtain the gradient of rotation angle at an arbitrary position. As we know in Mechanics of Materials, the structure that is subjected to torques or to bending moments can cause a three-dimensional rotation angle γ at an arbitrary position (x, y, z) of the structure, which is described in terms of the one-dimensional rotation angles γ i around the i-axis written as γ = γ x e x + γ y e y + γ z e z , (13) in which e i is a unit vector component around the i-axis. The one-dimensional rotation angles γ i are defined as in Eq. (1f).
Considering an arbitrary cube element abcdABCD with three lengths, i.e., dr of the line element AD, dr of aA and dr of AB (Fig. 2a), we suppose that one of end points position of vertical line element rs between the two parallel planes abBA and dcCD is (x, y, z) in plane abBA, and the other is (x + dx, y + dy, z + dz) in plane dcCD. Therefore, the vertical line element has the length vector dr and the direction cosine vector n respectively equal to where the operator represents the vector norm, and Similarly, the other two length vectors and the direction cosine vectors are given respectively for two parallel planes abcd and ABCD as dr = dx e x + dy e y + dz e z , n = υ x e x + υ y e y + υ z e z , (14d, e) and for aADd and bBCc as dr = dx e x + dy e y + dz e z , n = υ x e x + υ y e y + υ z e z . (14f, g) Subjecting the cube element abcdABCD to torques and to bending moments, therefore, we gain rotation angle components γ x (x, y, z), γ y (x, y, z) and γ z (x,y, z) of the plane abBA at point r, and γ x (x + dx, y + dy, z + dz), γ y (x + dx, y + dy, z + dz) and γ z (x + dx, y + dy, z + dz) of the plane cCDd at point s, where the line element rs is parallel to AD. Regarding the rotation angles of the plane abBA at point r as the references and on basis of the first order Taylor formula with the remain of Peano o( dr ), the components of the relative rotation angle vector between points r and s are The relative rotation angle vector dγ is projected onto the direction cosine vector n of the plane abBA, i.e., δ = dγ · n is relative twist angle between two parallel planes abBA and c C d d (Fig. 2b), and then is divided by the vector norm dr of the line element rs, and we get the gradient of twist angle ζ nn stated as where H(n) = dγ / dr called gradient of rotation angle vector is the relative rotation angle vector dγ per unit length along the direction n, and its components are given by Substituting Eqs (16b-e) and (14b) in Eq. (16a), one gives The relative rotation angle vector dγ is projected onto the direction n , i.e., ξ = dγ · n is relative bending angle between two parallel planes aADd and b B C c (Fig. 2c), and then is divided by the vector norm dr , and we get the gradient of bending angle α stated as where H(n ) is the gradient of rotation angle vector along the direction n and could gained by substituting n for n of H(n) in Eqs (16b-e). Similarly, we can also obtain the gradient of bending angle β between two parallel planes abcd and ABCD stated as β = H (n ) n .
(17b) Consequently, we can get the total gradient of bending angle ζ n n as ζ n n = α + β = H (n ) n + H (n ) n .
(17c) Substituting Eqs (14d-g), (16b-e) and (17a, b) in Eq. (17c), one gives Eqs (16f) and (17d) are respectively rewritten in tensor as The two tensor Eqs above are combined to a following tensor Eq. as whereη i j =η ji are symmetric curvature tensor. By using the Eqs (15) without the remain of Peano o( dr ), the rotation angle vector γ s at point s can be expressed in the block matrix form in terms of the rotation angle vector γ r at point r, symmetric curvature tensorη, skew-symmetric curvature tensorχ and length vector dr as γ s = γ r +η dr +χ dr, (19a) whereχ i j = −χ ji are skew-symmetric curvature tensor, which are given bỹ (19b-f) Let us define a gradient of rotation angle vector ω = −χ yz e x +χ xz e y −χ xy e z that consists of skew-symmetric curvature tensor χ i j , we can find that the orbital angle vector ω ×dr equals to the third term of right side of Eq. (19a), that is, Introducing another Cartesian coordinatexỹz whosez axis is parallel to the direction of rotation vector ω (Fig. 2d), we can rewrite ω and dr as ω = ω ez, dr = dxex + dỹeỹ + dzez, (21a, b) where ω is the vector norm of the gradient of rotation vector ω, and eĩ is a unit vector component around theĩ-axis.
Substituting Eqs. (21) in ω ×dr, we can obtain where For a small gradient of rotation angle ϖ corresponding to the gradient of rotation angle vector norm ω 1, the orbital angle vector norm ω ×dr approximatively equals to arc length bb as where arc length bb equals to the small gradient of rotation angle ϖ multiplied by length of segment sr , and sr is shadow of the line element sr inxỹ plane, and the following two equations can be obtained respectively Substituting Eqs (22b, d, e) in Eq. (22c), we can obtain the gradient of the rotation angle ϖ as = ω×dr When the line element sr rotates through the gradient of rigid angle ω , norm of line element sr is given by For the small rotation angle, i.e., ω 1, the above equation has an approximated relationship as Therefore, the vector ω ×dr suggests the end point s of the line element sr rotates through the gradient of rigid angle ω around the direction of rotation vector ω, where the rotation axis passes through the other end point r.
Substituting Eq. (20) in Eq. (19a), one gives γ s = γ r +η dr + ω×dr. (25) Above equation suggests that the rotation angle vector γ s of the plane cCDd at the point s consists of three parts: the rotation angle vector γ r of the plane abBA at the point r, rotation angle deformation vectorη dr and rigid rotation angle vector ω ×dr. Especially, symmetric curvature tensorη i j are two types of gradient of rotation angles: gradient of bending angleη i j (i = j) and gradient of twist anglẽ η ii (i = j); skew-symmetric curvature tensorχ i j are gradient of rotation angle for rigid body rotation.
In nonlinear elasticity, we also expect that the skew-symmetric curvature tensorχ i j will not contribute as a fundamental measure of deformation as the same as in the MCST and obtain the rotation angle vector γ s only in terms of the rotation angle vector γ r and the rotation angle deformation vectorη dr as γ s = γ r +η dr. (26) On basis of definition Eq. (16a) for the gradient of twist angle ζ nn , we can get the differential of twist angle γ n n between two arbitrary parallel plane elements along whose direction cosine n as dγ n n = dγ n = ζ nn dr .
Substituting Eqs (14c) and (18a) in above equation, we can obtain the differential of twist angle dγ n n along an arbitrary direction n as dγ n n = Thus, based on above equation, we also can obtain the differential of twist angle scalar γ i i along i-axis (i = x, y, z) as Supposing the gradient of twist angle scalar γ i,i constant, Eq. (27c) is rewritten by where ς i i are the weighting coefficients, and Integrating Eq. (27d) with respect to i from zero to L i under the condition that γ i i = 0 at i = 0, one gives where L i is structure length along i-axis.
In the same way, we also obtain bending angle scalar γ in which the gradient of bending angle scalar Considering a uniform circular cross-section shaft of length L i and of radius R i that was made of nonlinear elasticity materials and whose left end was fixed along i-axis (i = x, y, z), we exert a torque named T i parallel to i-axis at its the other end as shown in Fig.  3a. Due to deformation caused by the twisting couples, there is an angle of twist γ i i that is a one-dimensional twist angle parallel to i-axis.
Based on the assumption in Mechanics of Materials that the one-dimensional deformed circular shaft is made of separate slats parallel to the i-axis, its torsional deformation actually consists of shear deformation of these slats. We detach from the solid shaft an element of the hollow shaft of inner radius ρ i and of the outer ρ i + dρ i , and then consider the element of the hollow shaft made of element of separate slat (Fig. 3c), for example, the element of slat (Fig. 3b) that has element of line ab parallel to the axis of the shaft on the surface of shaft of inner radius ρ i and that has the element of area dA i on the cross-section of right end of the shaft (Fig. 3d).
As stated in Mechanics of Materials, the element of slat of the shaft will be subjected to shear stresses σ iγ and σ γ i , we could obtain the engineering shear strain ε γ i (i.e., the angle bab in Fig. 3b) corresponding to an angle of twist γ i i at the right end. For small values of shear strain ε γ i , we could both express the arc length bb as bb = L i ε γ i and as bb = Substituting Eq. (27f) in above equation, one gives The torques T i at the cross-section of right end of the shaft (Fig. 3d) is equal to the sum of the moments of shear forces σ γ i dA i rotating round the axis of the shaft, in other words, where ρ i denotes the vertical distance between the differential of shear force σ iγ dA i and the rotation axis of the shaft, and σ γ i is the shear stress on the element of area dA i as shown in Fig. 3d that is expressed as where dκ i is the element of central angle corresponding to the element of area dA i . Substituting Eqs (30) and (31b) into Eq. (31a) and integrating with respect to κ i from zero to 2π, and then with respect to ρ i from zero to R i , the torque T i is obtained as Now, let us consider the elementary work dλ i done by the torques T i as the shaft twists by one small amount dγ i i . The elementary work is expressed as Substituting Eqs (32) and (27f) in above equation and integrating with respect to i from zero to L i , the workλ i is obtained as Dividing the workλ i by the volume V = π(R i ) 2 L i of the undeformed circular shaft in Fig. 3a, we have the average strain-energy density i , i.e., the strain energy per unit volume, written as Taking the first derivative of Eq. (35) with respect to twist angle per unit length η ii , we could obtain torsion stress scalar m ii written as where and m α ii (α = L, N) represent linear and nonlinear components of the torsion stress scalar, respectively; l ii 1 and l ii 2 are material length scale parameters as stated in linear couple stress-strain constitutive relationships Eq. (1d).
Similarly, considering a uniform circular cross-section shaft that is subjected to bending moment parallel to i-axis, there are a bending angle γ j i and bending angle per unit length . We believe that through the derivation of the formulas like above, bending-moment stress scalar called the couple stress scalar m ij has the same expression as the torsion stress scalar m ii and is written as This is because torque consists of the sum of shear stress multiplied by its arm and bending moment consists of the sum of normal stress multiplied by its arm. What's more, for the nonlinear elasticity materials both shear stress and normal stress are nonlinear and can interconvert into each other.
On basis of Eqs (2a), (4b), (12), (27e), (28b) and (36), the strain energy density of an arbitrary volume element ϕ(ε ij ,η ij ) described in scalar is given by Where Setting the generalized Young's modulus coefficient E 2 = 0, the strain energy density Eq. (37a) will come back to the linear one. According to the Bernoulli-Euler beam theory [40], the displacement field can be written as

M ATHE M ATIC AL FOR MUL ATION
whereũ x (x, y, z, t ),ũ y (x, y, z, t ) andũ z (x, y, z, t ) are the displacements of an arbitrary space point (x, y, z) in the microbeams along the x, y and z directions, respectively; u(x, t) and w(x, t) represent the displacements of the mid-plane in the longitudinal and transverse directions, respectively; and t is time. Using Eq. (38) in Eqs (37i) and (37j), the nonzero components of η ij can be written as Similarly, substituting Eq. (38) into Eqs (37c) and (37e), the nonzero components of ε ij can be given as Taking into account von Kármán geometric nonlinearity due to the exact strain of large deformation, the nonlinear straindisplacement relationship [41,42] is obtained by By substituting Eqs (37b) and (37h) into Eq. (37a), the strain energy density can be obtained as Using Eqs (2b)-(2d) and (36f)-(36h), Eq. (41a) could be rewritten as follows: Consequently, the strain energy density equation (41b) based on the MCST in nonlinear elasticity has the additional nonlinear normal stress σ N xx and couple stresses m N yx compared with Eq. (1a) based on the MCST in linear elasticity. Next, by substituting Eqs (39) and (40b) into Eq. (41b), the strain energy of the isotropic nonlinear elasticity material microbeam occupying region can be derived as where the normal resultant force N xx , bending moment M xx and couple moment Y yx , and the nonlinear generalized normal resultant force N N xx , bending moment M N xx and couple moment Y N yx , which are caused by the nonlinear elasticity, can be defined as follows: The corresponding coefficients in Eqs (42b)-(42g) are given by For the simplification of Eq. (42a), one obtains where On the other hand, the kinetic energy of the system can be expressed as where The governing equations for the free vibration of microbeams can be derived from Hamilton's principle [9,43]: where δ is the first variation operator and F is the Lagrange function of the free vibration system without the external forces, which is given by (45c-g) On the basis of the derivation rule of compound function, the first variation of Lagrange function can be rewritten as Substituting Eqs (43a), (44a) and (45b) into Eq. (46) and setting the coefficients of δw, δu and δw , x to zero yield the nonlinear partial differential governing equations of microbeams: the transverse equation the longitudinal equation and the three kinds of boundary ends: C end We obtain the Lagrange function F, which could be expressed in derivatives of the displacement field with respect to x and t (because N*, M* and Y* are functions of derivatives of the displacement field and F is compound function of N*, M*, Y* and derivatives of the displacement field), and then take the variation of the Lagrange function F with respect to derivatives of the displacement field. Thus, there are many derivatives of N*, M* and Y* with respect to derivatives of the displacement field (such as ∂N * /∂w ,x ) in the abovementioned equations. Meanwhile, there are no derivatives of the normal resultant force, bending moment and couple moment with respect to derivatives of the displacement field as obtained byŞimşek et al. [42,44]. This is because they take the variation of the strain energy density with respect to the strain directly and then take the variation of the Lagrange function F with respect to derivatives of the displacement field. The two procedures could obtain the same results even though a difference in expression exists.

SOLU TION TECHNIQUE
To gain the dimensionless equations of motion, the dimensionless parameters are defined as follows: where L 0 indicates that the length L of microbeams is constant. On the basis of the discretization rules of the DQM [27], the dimensionless displacement field components W, U and their nth derivatives with respect to X can be approximated as where N is the total number of nodes distributed along the x direction, γ m (X) are the Lagrange interpolation polynomials and C (n) im are the weighting coefficients. The sampling points are generated by using the following cosine pattern: Substituting Eqs (51) and (52) into Eqs (47)-(50), the discretized ordinary differential equations governing the nonlinear vibration of microbeams are written in the block matrix form as ⎡   Equations (54) and (55) are rewritten in more simplified matrix form for obtaining eigenvalue equations easily as (57) Inserting Eq. (56) into Eq. (57), one obtains  Next, the dynamic displacement vector {D} b can be expanded in the form of Then, substituting Eq. (59a) into Eq. (58a) yields the nonlinear eigenvalue equations as shown below: The above nonlinear eigenvalue equation (60) can be solved through a direct iterative process as follows, applied in [14,15,20,21]: Step 1: By neglecting nonlinear terms in the stiffness matrix [K] of Eq. (60), we could obtain the linear stiffness matrix [K] L and the corresponding linear eigenvalue problem is solved.
Step 2: The linear eigenvectors obtained in Step 1 are appropriately scaled up such that the maximum transverse displacement is equal to a given vibration amplitude. Then, inserting the scaled normalized linear eigenvectors into the stiffness matrix [K], we could use the scaled normalized linear eigenvectors to get the nonlinear stiffness matrix [K] NL . The nonlinear eigenvalues and eigenvectors are obtained from the updated eigensystem (60).
Step 3: The eigenvectors are scaled up again, and Step 2 is repeated until the frequency values from the two subsequent iterations i and i + 1 satisfy the prescribed convergence criteria as where ω k is the frequency at iteration k(k = i, i + 1) and ε 0 is a small value number, which is set to be 10 −4 in this paper.

NUMERICAL RESULTS
Numerical results of the nonlinear free vibration of microbeams made from nonlinear elasticity material Al-1% Si are presented in Tables 1-4 and Figs 5-13. In previous studies [36,38], two generalized nonlinear Young's modulus coefficients of Eq. (3a) of Al-1% Si were given as E 1 = 65 GPa and E 2 = 6 × 10 5 GPa, and the linear elastic modulus was E = 69 GPa. However, the generalized Young's modulus coefficient E 1 of Eq. (2a) should be equal to E of Eq. (1g) when ε xx = 0. To obtain the relationship E 1 = E, we use the third-order nonlinear algebraic equation of Eq. (2a) again to refit the nonlinear normal stress-strain constitutive relationship   or linear (L) classical stress, and the third subscript ς represents nonlinear (N) or linear (L) couple stresses. For example, ω 1 LLN is the fundamental nonlinear dimensionless frequency within the framework of the geometric linearity, and the linear classical and nonlinear couple stress-strain constitutive relationships; W 1 NNN is the fundamental nonlinear dimensionless mode shape within the framework of the geometric nonlinearity, and the nonlinear classical and couple stress-strain constitutive relationships. Table 1    To validate the present study, two direct comparisons with theoretical predictions of the linear free vibration responses in Table 2 and numerical results of the nonlinear free vibration responses in Table 3 are made, respectively. A size effect exists in the MCST when l = 0. Otherwise, no size effect exists in the classical beam theory. Table 2 lists the first three linear natural frequencies ϕ j LLL (j = 1, 2, 3) of microbeams with four groups of boundary conditions on the basis of geometric linearity, and Table 3 displays the fundamental nonlinear frequency ratio ϕ 1 NLL /ϕ 1 LLL of microbeams with three groups of boundary conditions on the basis of geometric  Ke et al. [14,15] found that employment of von Kármán geometric nonlinearity for nonlinear free vibration exhibits a typical "hard spring" behavior, i.e. taking into account von Kármán geometric nonlinearity increases vibrational frequencies as the vibration amplitude increases. Therefore, to eliminate the disturbance of von Kármán geometric nonlinearity in Eq. (40b), the linear strain-displacement relationship, i.e. Eq. (40b) without von Kármán geometric nonlinearity term (∂w/∂x) 2 /2, is adopted to research the effect of nonlinear elasticity properties on free vibration behaviors of microbeams. Figures 5-7 plot the comparisons between the first three nonlinear dimensionless frequencies ω j LLN (j = 1, 2, 3) and linear dimensionless frequencies ω j LLL with the H-H, C-C, C-H and C-F ends under different dimensionless amplitudes W 1 max (H/l = 1, L/H = 100). At a given positive amplitude W 1 max , the nonlinear frequencies ω j LLN with nonlinear couple stress are smaller than linear frequencies ω j LLL with linear couple stress. In addition, microbeams with nonlinear couple stress-strain constitutive relationships exhibit a typical "soft spring" behavior, i.e. the nonlinear frequencies ω j LLN with nonlinear couple stress decrease gradually as the vibration amplitude W 1 max increases.  [36][37][38] obtained the same results by studying the forced vibrations of microbeams with nonlinear classical stress-strain constitutive relationships. In addition, we find microbeams with nonlinear classical stress-strain constitutive relationships also exhibit a typical "soft spring" behavior, i.e. the nonlinear frequencies ω j LNL including nonlinear classical stress decrease gradually as the vibration amplitude W 1 max increases. Figures 11-13 show the comparisons between the first three nonlinear dimensionless frequencies ω  Table 4, where the percentage values in parentheses represent (ω It is clear that all the relative error values are negative, which implies that taking into account both nonlinear classical stress and couple stress decreases the frequencies. At the given vibration amplitudeW 1 max = 1, the C-C microbeams have the highest absolute values of relative error of the first three frequencies (-4.57%, -1.11% and -1.46%) and the C-F microbeams have the lowest values (-0.092%, -0.028% and -0.029%). These results show that taking into account both nonlinear classical stress and couple stress has a significant effect in softening stiffness.

CONCLUSIONS
In this study, the nonlinear couple stress-strain constitutive equations are obtained on the basis of previous nonlinear stress-strain constitutive relationships of nonlinear elasticity materials. Hamilton's principle is utilized to derive higher-order nonlinear governing equations and boundary conditions within the framework of the updated MCST, von Kármán geometric nonlinearity and Bernoulli-Euler beam theory. The nonlinear partial differential governing equations are discretized into a set of nonlinear ordinary differential equations using the DQM. Then, an iterative algorithm is employed to solve these nonlinear ordinary differential equations for obtaining the nonlinear vibration frequencies of the microbeams with four types of boundary conditions. The influences of nonlinear elasticity properties on the nonlinear free vibration characteristics of the microbeams are discussed. The following conclusions are obtained: (1) the use of nonlinear couple constitutive relationships for nonlinear elasticity materials has a significant softening effect on the stiffness, i.e. employing nonlinear couple stress-strain constitutive relationships leads to lower vibrational frequencies than the previous linear couple stress-strain constitutive relationships; and (2) the nonlinear elasticity materials exhibit a typical "soft spring" behavior, i.e. their stiffness softens as the vibration amplitude increases.
Expressions of block matrices in Eq. (57) are written for governing equations as follows: