On Strain Gradient Theory and Its Application in Bending of Beam

: The general strain gradient theory of Mindlin is re-visited on the basis of a new set of higher-order metrics, which includes dilatation gradient, deviatoric stretch gradient, symmetric rotation gradient and curvature. A strain gradient bending theory for plane-strain beams is proposed based on the present strain gradient theory. The stress resultants are re-deﬁned and the corresponding equilibrium equations and boundary conditions are derived for beams. The semi-inverse solution for a pure bending beam is obtained and the inﬂuence of the Poisson’s effect and strain gradient components on bending rigidity is investigated. As a contrast, the solution of the Bernoulli–Euler beam is also presented. The results demonstrate that when Poisson’s effect is ignored, the result of the plane-strain beam is consistent with that of the Bernoulli–Euler beam in the couple stress theory. While for the strain gradient theory, the bending rigidity of a plane-strain beam ignoring the Poisson’s effect is smaller than that of the Bernoulli–Euler beam due to the inﬂuence of the dilatation gradient and the deviatoric stretch gradient along the thickness direction of the beam. In addition, the inﬂuence of a strain gradient along the length direction on a bending rigidity is negligible.


Introduction
A micro-beam is a typical component in micro-electro-mechanical systems (MEMSs).Accurate description of mechanical behavior of a beam is necessary for the design and performance prediction of MEMSs [1][2][3].Although the traditional beam model has matured, the mechanical behavior of beams at the micro-scale is obviously different from that at the macro-scale, which cannot be captured by the traditional models [4][5][6].Many experiments have reported that the mechanical properties of a micro-beam show obvious size effects [7][8][9][10].McFarland and Colton observed that the bending stiffness of polypropylene cantilever beams with thicknesses of 15 and 30 µm is more than four times the value predicted by traditional theory [9].Lam et al. found that the dimensionless bending stiffness of the epoxy polymeric beam with a thickness of 20 µm was approximately 2.3 times higher than that of 115 µm [7].For the purpose of capturing the size-dependent mechanical properties, the strain gradient theory has been presented.
In strain gradient theory, the strain gradient is regarded as an additional (higher-order) deformation metric in addition to the traditional strain.Among the strain gradient theories, the couple stress theory is the earliest strain gradient theory that uses the rotation gradient as an additional deformation metric [11].In fact, the stain gradient includes dilatation gradient, rotation gradient and deviatoric stretch gradient [12].The rotation gradient is only the anti-symmetric part of the strain gradient.As with the development of strain gradient theory, the general strain gradient theory, including all strain gradient components, is proposed by Mindlin and Eshel by introducing five higher-order material constants [13].Then, the general strain gradient theory of Mindlin has been re-formulated by Lam et al. on the basis of the equilibrium of moments of couples [7], and by Zhou et  independent strain gradient components [14], respectively.The two re-formulated strain gradient theories include three higher-order material constants.A simple strain gradient theory with only one higher-order material constant has been presented by Aifantis [15].The general form of couple stress theory includes two higher-order material constants, which has been further simplified to include only one higher-order material constant by Yang et al. [16] and Hadjesfandiari and Dargush [17], respectively, by choosing only the symmetric or anti-symmetric parts as the higher-order deformation metrics.However, in our opinion, although these re-formulated strain gradient theories with three or fewer higher-order constants can also capture the size-effect phenomenon, these theories ignore the contribution of some strain gradient components and underestimate the stiffening effect of the strain gradient.In fact, Mindlin's strain gradient theory is a general strain gradient theory and should be used to develop micro-component models.
According to the strain gradient theory, many strain gradient beam models have been established to capture the size-dependent mechanical behaviors [18][19][20][21][22][23].Based on Mindlin's strain gradient theory, a linear model of Bernoulli-Euler beams has been presented by Niiranen et al. [24], and a geometrically non-linear model was further developed by Tran and Niiranen [25] by adopting the von Kármán strain assumption.Based on the modified couple stress theory, Vo et al. presented the models of planar [26] and spatial [27] arbitrarily curved micro-beams and Li et al. investigated the bending and free vibration of bi-directional functionally graded graphene nano-platelets-reinforced composite microbeams [28].However, these existing models are contradictory [29].The non-zero strain of Bernoulli-Euler beams is an axial strain, and its non-zero strain gradient includes the gradients of axial strain along the thickness direction and along the length direction.Currently, there are three main types of models available, depending on the stain gradient components used.The first category considers the strain gradient along the thickness direction and the length direction at the same time [30][31][32][33].The second type considers only the strain gradient along the length direction and does not consider the strain gradient along the thickness direction [34].The third type considers the strain gradient along the thickness direction rather than the length direction, since the strain gradient along the length direction can be neglected compared with that along the thickness direction for Euler beams [35].The third type may be seen as an approximation of the first type.However, the second and third types are polar opposites, and one of them must be wrong.Thus, it is important to figure out which type is correct.
Lurie and Solyaev [29] checked the beam models of the first and second types and believed that only the second type (uniaxial stress state) could guarantee the satisfaction of boundary conditions on the top and bottom surfaces of the beams.However, in our opinion, the satisfaction of boundary conditions on the top and bottom surfaces of the beams depends on the definition of moment, which includes the higher-order stress.The first type can satisfy the boundary conditions on the top and bottom surfaces of the beams by defining a different moment from that of Lurie and Solyaev.
The purpose of this work is to establish a correct beam model and analyze the contribution of the strain gradient components to bending stiffness.This paper is organized as follows: Mindlin's general strain gradient elasticity theory is re-formulated by introducing a new set of strain gradient components so as to investigate the effects of each strain gradient component conveniently in Section 2.Then, the bending theory for plane-strain beams is presented in Section 3, and the semi-inverse solution for pure bending beams is developed in Section 4. Section 5 investigates the bending of Bernoulli-Euler beams and compares the results of plane-strain beams and Bernoulli-Euler beams.Finally, conclusions are summarized in Section 6.

A Summary of Strain Gradient Theory
In strain gradient theory, not only the traditional strain ε ij but also its gradients η ijk are considered in the internal energy density U, written as: with where u i denotes the displacement vector and a comma represents the differentiation with respect to the coordinates.The work-conjugated Cauchy stress σ ij and higher-order stress τ ijk with respect to strain and strain gradient, respectively, are given by: Obviously, the stress tensor σ ij is a symmetric tensor and the higher-order stress tensor τ ijk is symmetric with respect to subscripts j and k.
For a solid of volume V with surface boundary S and sharp edge C, the virtual work principle of the strain gradient theory is given as [14]: Thus, the equilibrium equation can be obtained by applying the variation of internal energy, as shown in Appendix A in detail, given as: and the boundary conditions on surface boundary S and sharp edge C are: and where D = n i ∂ i is the normal gradient operator with n i representing a unit vector normal to the boundary surface S, D i = (δ ik − n i n k )∂ k represents the surface gradient operator, the square brackets stands for the difference between the values of the enclosed quantity on the two sides of the edge and k i denotes the outer co-normal vector.
For linear elastic isotropic materials, the internal energy density of Equation ( 1) is generally written as: where ε ′ ij = ε ij − 1 3 δ ij ε nn denotes the deviatoric strain tensor.Then, the constitutive equations in terms of Equation (3) are given by: Coatings 2022, 12, 1304 4 of 19 in which k and µ are the bulk and shear moduli, δ ij denotes the Kronecker delta and g i (i = 1, 2, 3, 4, 5) is the introduced higher-order constant.

Development of Strain Gradient Components, and Equilibrium and Boundary Conditions
Strain gradient has been decomposed into its orthogonal components to re-formulate the strain gradient theory [14].In fact, strain gradient includes dilatation gradient γ i , deviatoric stretch gradient η (1) ijk and rotation gradient χ ij .These strain gradient components are defined as: with e ilk denoting the alternating tensor.The rotation gradient can be further decomposed into its symmetric and anti-symmetric parts, given as: Additionally, the anti-symmetric rotation gradient tensor χ a ij can be represented by its dual vector κ i which is called as the curvature vector [17], written as: Then, the strain gradient can be finally expressed as the sum of dilatation gradient γ i , deviatoric stretch gradient η ijk , symmetric rotation gradient χ s ij and curvature κ i , given as [12]: with the coefficients being: Thus, the internal energy density becomes a function of strain and strain gradient components, written as ijk , χ s ij , κ i ).The corresponding work-conjugated higher-order stress components p i , τ ijk , χ s ij , κ i , respectively, are given by: The virtual work density can be rewritten in terms of the new strain metrics as: ijk δη (1) Coatings 2022, 12, 1304 5 of 19 From Equation ( 18), it can be concluded that the higher-order stress τ ijk is associated with the higher-order stress components p i , τ ijk , m s ij , µ i by: Finally, applying the principle of virtual work yields the equilibrium equation in V as: and the boundary conditions on S as: and along the sharp edge C as: For linear elastic isotropic materials, the internal energy density in Equation ( 9) can be re-formulated by substituting Equation (15) into Equation ( 9), written as: ijk η (1) where a i (i = 1, 2, 3, 4, 5) is a higher-order material constant in another form which is related to g i by: By re-writing the higher-order material constants as: with l i (i = 0, 1, 2, 3, 4) representing five material length scale parameters, the higher-order constitutive relations can be obtained from Equations ( 17) and (25) as: The present strain gradient theory is another expression of Mindlin's strain gradient theory, which is a general theory with five higher-order material constants.When the curvature vector κ i is excluded (by letting l 3 = l 4 = 0), the present formulation of strain Coatings 2022, 12, 1304 6 of 19 gradient theory will reduce to that of Lam et al.'s [7].When the material length scale parameters are constrained as: the internal energy density of the present strain gradient theory will reduce to which is the internal energy density of Zhou et al. with χ ′ ij = e ilk ε ′ jk,l denoting the deviatoric rotation gradient [14].

Couple Stress Theory
In couple stress theory, the rotation gradient rather than the strain gradient is included in the internal energy density.That is, the strain gradient theory will reduce to the couple stress theory when the dilatation gradient γ i and deviatoric stretch gradient η (1) ijk are excluded.Thus, when the higher-order stress components p i and τ (1) ijk are excluded, the equilibrium equation and boundary conditions of the strain gradient theory in Equations ( 20)-( 24) can reduce to those of the couple stress theory.The equilibrium equation is: and the boundary conditions on S are: and along the sharp edge C is: In addition, the anti-symmetric couple stress component m a ij , which is work-conjugated with the anti-symmetric rotation gradient component χ a ij , is commonly used in many articles rather than the present higher-order stress component µ i .In fact, according to the virtual work density of couple stress theory The anti-symmetric couple stress component m a ij is associated with the higher-order stress component µ i by m a ij = 1 2 e kji µ k .Additionally, in turn, the higher-order stress component µ i is related to the anti-symmetric couple stress component m a ij through µ i = e ijk m a kj .Then, from Equations ( 31)- (34), the equilibrium equation and boundary conditions of couple stress theory can be re-formulated in the form of couple stress components m s ij and m a ij ,which are shown in Appendix B in detail.
For linear elastic isotropic materials, the internal energy density of the couple stress theory can be reduced from Equation (25), written as: Consider the material length scale parameters in Equation ( 27), the couple stress components can be given as:

Bending Theory for Plane-Strain Beams
Consider the same plane-strain beam problem as Lam et al. [7].The length of a beam is L and the thickness is h.The thickness of a beam is much smaller than its length.The Cartesian coordinate system is adopted, and the X-axis coincides with the center line of the beam along the length direction, and the Z-axis coincides with the thickness direction.By ignoring the body forces and only considering the normal tractions q ±h/2 in units of N/m 2 applied on the top or bottom surfaces of the beam, the governing equations along the length and thickness directions can be obtained from Equation (20), given by: The boundary conditions on the top and bottom surfaces of the beam, z = ±h/2, from Equations ( 21)-( 23) are: 333,3 = q ±h/2 (40) On the surfaces normal to the X-axis, x = 0 or L, the boundary conditions obtained from Equations ( 21)-( 23) are: ,1 113 + (47) The boundary conditions at the edges of the top and bottom surfaces of the beam from Equation ( 24) are: The stress resultants can be re-defined due to the presence of higher-order stress components as: where N, Q, M, N h and M h denote the axial stress resultant, shear stress resultant, moment, higher-order stress resultant and higher-order moment per unit width in the beam, respectively.According to the re-defined stress resultants in Equation ( 49) and the governing equations and boundary conditions in Equations ( 38)-( 48), the same formulations from that of Lam et al. [7] can be obtained as: where dQ dx + q = 0 (51) with q = q h/2 + q −h/2 .Substituting Equation ( 52) into (51) yields: The stress resultants at the boundary are: where N, Q, M, N h and M h represent the prescribed axial force, shear force, moment, higher-order axial force, higher-order moment per unit width in the beam, respectively.

Semi-Inverse Solution for Pure Bending Beam
Consider a pure bending beam subjected to a moment M 0 acting on its two ends.The governing equations and boundary conditions are listed as Equations ( 38)-(48) with normal tractions q ±h/2 being zero.The solution for the displacement field is suggested in the following form [29]: Coatings 2022, 12, 1304 9 of 19 in which w is an unknown function with respect to coordinate z.C 1 , C 2 , u 0 and v 0 are four undetermined constants.According to the assumption of displacement, the non-zero strain components are: The non-vanishing components of strain gradient are calculated as: According to the constitutive relations for plane-strain problems, the non-zero Cauchy stress is given by: in which E is the Young's modulus and ν is the Poisson's ratio.The work-conjugated higher-order stress components can be obtained from Equation (28), given as: According to the present bending theory of plane-strain beams, the governing equations and boundary conditions of the pure bending beam subjected to a moment M 0 acting on its two ends can be specified by substituting the stress components from Equation (62) and higherorder stress components from Equation (63) into Equations ( 38)-(48) along the length direction of Equation ( 38) and the boundary conditions of Equations ( 40), ( 42), ( 45) and ( 47) can be satisfied automatically.The current pure bending solution depends on the governing equations, which are Equation (39) and boundary conditions Equations ( 41), ( 43), ( 44), ( 46) and (48).
Substituting the stress components from Equation (62) and higher-order stress components from Equation (63) into Equation (69) yields the governing equation of the current pure bending beam, that is: with the coefficient d 1 defined as The solution of Equation ( 64) is given by: where C 3 is an undetermined coefficient and the parameter λ is defined as: It should be noted here that the constant term is not included in Equation (65) because it has been included in v 0 in the displacement solution of Equation (59).
Moreover, the assumption of displacement in Equation (65) satisfies the boundary condition in Equation (41) automatically.The coefficients C 1 and C 3 in Equation (65) can be determined from the boundary conditions Equations ( 43), (44), ( 46) and (48).Based on the boundary condition Equation ( 43), it can be obtained that  56), actually, the applied moment M 0 at the two ends of beam equals Then, combining Equations ( 44), ( 46), ( 48) and ( 68), the coefficient C 1 can be solved as: with the coefficients d 3 and d 4 denoting Furthermore, the coefficients C 2 , u 0 and v 0 can be obtained according to the displacement boundary conditions.For a simply supported beam, the displacement boundary conditions are: Substituting Equations ( 59) and ( 65) into (71) yields Thus, the displacement solution is finally given by: in which the coefficient C 1 is given as Equation (69).Especially, for the axis of the beam (z = 0), its displacement solution is: which is similar to that of traditional theory except for the bending rigidity.From Equation (69), it can be seen that the bending rigidity of the pure bending beam based on the strain gradient theory includes the classical bending rigidity Eh 3 12(1−v 2 ) and the higher-order bending rigidity . By using the classical bending rigidity to normalize the total bending rigidity of the strain gradient theory (the denominator in Equation ( 69)), the dimensionless bending rigidity K can be obtained as: with the equivalent material length scale parameter l denoted as: Coatings 2022, 12, 1304 11 of 19 When the Poisson's effect is ignored by ν = 0, the dimensionless bending rigidity K reduces to in which the parameter λ is reduced from Equation (66) by letting v = 0.
Similarly, for a couple stress plane-strain beam subjected to a moment M 0 acting on its two ends, the displacement solution based on the couple stress theory can be obtained as: where the coefficient is: Then, the dimensionless bending rigidity for the couple stress theory is given as: with the equivalent couple stress length scale parameter l c denoted as l c = l 2 2 + l 2 3 .When the Poisson's effect is ignored, the dimensionless bending rigidity K c reduces to In order to illustrate the solution of the present plane-strain beam, a polyvinylidene difluoride (PVDF) beam with a unit thickness is considered here, in which the Young's modulus E is equal to 3.7 GPa, the material length scale parameters are assumed to be the same, the length of the beam is 20 times its thickness (L = 20 h), and the applied moment M 0 equals 10 µN•µm.
Based on Equations ( 75) and (81), the dimensionless bending rigidity is shown in Figure 1.It can be seen from Figure 1 that the dimensionless bending rigidity shows an obvious size effect.The dimensionless bending rigidity of beams decreases with an increase in thickness, and the bending rigidities based on the strain gradient theory, couple stress theory and classical theory tend to be consistent as the thickness increases.The effect of the strain gradient is a kind of stiffening effect.The bending rigidity based on the couple stress theory is larger than that of classical theory due to the consideration of a rotation gradient.Additionally, the bending rigidity based on the strain gradient theory is larger than that of the couple stress theory, since the dilatation gradient and deviatoric stretch gradient are also considered in addition to the rotation gradient.In addition, the dimensionless bending rigidity of a plane-strain beam based on higher-order theories increases with the decrease of the Poisson's ratio.
Figure 2 shows the displacement u 3 of the axis (Equation (74)) along the length of the beam.It can be found that the material length scale parameter l 0 , which characterizes the dilatation gradient, exhibits a dominant stiffening effect, while the material length scale parameters l 1 , l 2 and l 3 characterize the deviatoric stretch gradient, symmetric rotation gradient and curvature, respectively, exhibiting a secondary and decreasing stiffening effect.In particular, the material length scale parameter l 4 exhibits a softening effect.It should be noted here that the overall effect of the strain gradient is a stiffening effect, and the influence of l 4 is weak.-ε Moreover, in order to reveal the deformation of a plane-strain beam along the thickness direction, the strain ε 33 can be obtained by combining Equations ( 2) and (73) from the strain gradient theory, which is given by: While for couple stress theory, the strain ε 33 is obtained by: By defining the normalized strain ε 33 /C 1 for the strain gradient theory and ε c 33 /C c 1 for the couple stress theory, the normalized strain of a plane-strain beam along its thickness direction is shown in Figure 3. Figure 3 reveals that the gradient of a strain ε 33 along the thickness direction decreases as the Poisson's ratio decreases.When the Poisson's effect is neglected, the strain ε 33 is equal to zero for the couple stress theory and non-zero for the strain gradient theory.Equation (83) demonstrates that the strain ε 33 of a plane-strain beam based on the strain gradient theory includes the linear term and the non-linear term.The non-linear term is related to the parameter λ, which characterizes the dilatation gradient and the deviatoric stretch gradient along the thickness direction.When the dilatation gradient and the deviatoric stretch gradient along the thickness direction are neglected by letting λ = 0, the non-linear term will disappear.

Strain Gradient Bernoulli-Euler Beam
Consider a Bernoulli-Euler beam with a unit thickness, the displacement components are set as: According to the definitions of strain and strain gradient in Equation ( 2), the non-zero components are: Hence, the non-vanishing strain gradient components are: From Equation (49), the defined bending moment and higher-order moment are given as where the coefficients d 1 and d 3 have been defined in Section 4. The bending governing equation of the beam is given as Equation ( 53) and the boundary conditions of stress resultants are listed as Equations ( 55), ( 56) and ( 58).
(1) Simply supported beam For a simply supported beam subjected to a moment M 0 acting on the two ends of beam, the governing equation is: and the boundary conditions on the two ends of beam are: Combining Equations ( 89)-( 91), the deflection of pure bending beam is solved as: where the coefficient is given by: By using the classical bending rigidity Eh 3 /12 to normalize the bending rigidity based on the strain gradient theory (the denominator in Equation ( 93)), the dimensionless bending rigidity K B can be obtained as: where the equivalent material length scale parameter l is defined as Equation (76).In addition, when the dilatation gradient γ i and deviatoric stretch gradient η (1) ijk are excluded, the strain gradient pure bending beam model will reduce to that of the couple stress theory.Additionally, the dimensionless bending rigidity will be the same as Equation (82).
The displacement distribution of the pure bending beam along its axis is shown in Figure 4. Figure 4 illustrates once again that the bending rigidity predicted by the strain gradient theory is greater than that predicted by the couple stress theory, and the bending rigidity decreases with the increase of the Poisson's ratio.Moreover, for the couple stress theory, the displacement of the Bernoulli-Euler beam is consistent with that of the planestrain beam when the Poisson's effect is ignored.While for the strain gradient theory, the bending rigidity of the plane-strain beam without the Poisson's effect is smaller than that of the Bernoulli-Euler beam due to the influence of the dilatation gradient and the deviatoric stretch gradient along the thickness direction, which can also be explained by Equations (77) and (94).(2) Cantilever beam For a cantilever beam subjected to a concentrated force F 0 acting on its free end, the governing equation is the same as Equation ( 90) and the boundary conditions on the two ends of beam are given by: w = w ,1 = 0, M h = 0 at x = 0 (95)

   
Combining Equations (90), ( 95) and (96), the deflection is solved as: where λ and K B are defined as Equations ( 78) and (94), respectively, and the coefficients B i (i = 1,2,3,4,5) are: In addition, many articles only consider the gradient of the axial strain along the thickness direction since the strain gradient along the length direction can be neglected compared with that along the thickness direction for the Euler beam.For a cantilever beam subjected to a concentrated force at its free end, when the strain gradient along the length direction is neglected, the deflection solution of Equation (97) will reduce to w = B 3 x 2 + B 4 x 3 (99) From Equations ( 97) and (99), the deflection of the cantilever beam is shown in Figure 5.It can be seen that the influence of the strain gradient along the length direction η 111 on the bending rigidity is negligible.The main stiffening effect comes from the strain gradient along the thickness direction.In addition, compared with the couple stress model, the contribution of the dilatation gradient and deviatoric stretch gradient is non-negligible.

Conclusions
In this paper, a new set of higher-order deformation metrics including dilatation gradient, deviatoric stretch gradient, symmetric rotation gradient and curvature is introduced to re-formulate the strain gradient theory for the sake of investigating the effect of each strain gradient components conveniently.All the constitutive equations, equilibrium equations and boundary conditions in the form of components are obtained.The present strain gradient theory with five higher-order constants is another form of Mindlin's strain gradient theory, which is a general theory.The present theory can be directly simplified to the couple stress theory by eliminating the dilatation gradient and deviatoric stretch gradient.
Based on the present strain gradient theory, the bending theory of plane-strain beam is proposed, in which the stress resultants are re-defined due to the presence of higher-order stress and the corresponding governing equations and boundary conditions of stress resultants are developed.The semi-inverse solution of a pure bending beam subjected to moments at its two ends is obtained.The results demonstrate that the material length scale parameter l 0 , which characterizes the dilatation gradient, exhibits a dominant stiffening effect, while the material length scale parameters l 1 , l 2 and l 3 characterizing the deviatoric stretch gradient, symmetric rotation gradient and curvature, respectively, exhibit a secondary and decreasing stiffening effect.The dimensionless bending rigidity increases with the decrease in Poisson's ratio.
The Bernoulli-Euler beam model is presented based on the present strain gradient theory and couple stress theory.For the couple stress theory, the displacement of the Bernoulli-Euler beam is consistent with that of a plane-strain beam without the Poisson's effect.While for the strain gradient theory, the bending rigidity of a plane-strain beam without considering the Poisson's effect is smaller than that of a Bernoulli-Euler beam, which is due to the influence of the dilatation gradient and the deviatoric stretch gradient along the thickness direction.In addition, it can be concluded that the influence of the strain gradient along the length direction on bending rigidity is negligible.

Figure 2 .
Figure 2. The displacement u 3 of the axis along the length direction of the beam.

Figure 3 .
Figure 3. Normalized strain of a plane-strain beam along the thickness direction.

Figure 4 .
Figure 4.The displacement of beam along its axis.
al. according to the