# Numerical Research on a Three-Dimensional Solid Element Based on Generalized Elasticity Theory

Authors: Zhen Qiao; Zhanfang Liu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

A finite element equation is established based on generalized elasticity theory by applying a virtual work principle. Then, a penalty function term is added to the virtual work equation by imposing rotation and displacement as independent variables. An 8-node element with full integration, an 8-node element with reduced integration, and a 20-node element with full integration are constructed using difference integration schemes and shape functions. The influences of structural degrees of freedom and the penalty parameter on convergence are analyzed via the three elements. It is shown that the 8-node element with reduction integration and the 20-node element with full integration are convergent, whereas the 8-node element with full integration is divergent. The scale effects of a slender beam, a short beam, a thin plate, and a medium-thick plate are numerically analyzed. Lastly, the scale effects of the frequencies that correspond to the bending mode, torsion mode, and tension-compression mode for a pretwisted plate are studied. It is found that the frequencies that correspond to the bending mode and torsion mode exert a scale effect, whereas the frequency that corresponds to the tension-compression mode does not. The essence of the scale effect is that the rotational deformation of the microstructure is amplified.


---

## Body

## 1. Introduction

Microelectromechanical systems (MEMS) are primarily used in microsensors, microactuators, microresonators, and microscopes, in which beams [1–6] and plates [7–12] are the key models for analyzing MEMS. The mechanical properties of microstructures have scale effects, which had been observed in many experiments [13, 14]. Thus, this effect should be considered in the design and fabrication of MEMS. However, classical elasticity mechanics cannot measure the scale effect. Accordingly, several theories such as couple stress theory, strain gradient theory, and micropolar theory have been proposed to interpret the scale effect. Among these theories, couple stress theory has elicited considerable attention because of its simple constitutive model and few material parameters. Different from the classical elasticity theory (CET), the material particle in couple stress theory has six degrees of freedom (DOFs), i.e., three displacements plus three microrotations associated with displacement, and the Cauchy stress is asymmetrical. In addition, the constitutive model of the couple stress theory is controversial. Mindlin and Tiersten have proposed a linear constitutive model with four material parameters (two Lamé parameters and two additional material parameters) [15]. Subsequently, Mindlin developed a simpler constitutive model with only one additional material parameter [16]. However, these two constitutive models do not address the issue of the symmetry of the couple stress tensor. Yang et al. [17] argued that the couple stress tensor is symmetric and proposed a new constitutive model with one additional material parameter. By contrast, Hadjesfandiari and Dargush [18] believed that the couple stress tensor is antisymmetric and presented a new constitutive model with one additional material parameter. Huang and Huang [19] and Liu et al. [20] asserted that the couple stress tensor is asymmetric, and the latter renamed Mindlin’s couple stress theory [16] as generalized elasticity theory (GET).The displacement gradient that characterizes infinitesimal deformation can be divided into a symmetric strain tensor, which describes translational deformation, and an antisymmetric rotational tensor, which describes rotational deformation. The rotational tensor does not measure the degree of rotational deformation. Therefore, a curvature tensor is introduced to measure the degree of such deformations. Rotational deformation plays a vital role in microelastic bodies but exerts a minimal effect on macroelastic bodies. Therefore, GET, which retains rotational deformation, can measure the scale effect, whereas CET, which disregards rotational deformation, cannot measure such an effect. In GET, the relation between stress and strain describes translational deformation, whereas the relation between couple stress and curvature measures rotational deformation.The analytical solutions for many generalized elasticity problems are difficult to obtain. Therefore, numerical methods may be a better solution. Among numerical methods, the finite element method (FEM) is relatively mature and widely used. Accordingly, it can be adopted as a solution for generalized elasticity problems. Some elements that can be used in certain simplified models or specific situations have been mentioned in the literature in simulating the scale effect. These elements include plane elements [21, 22], hybrid elements [23, 24], Euler–Bernoulli beam element [25], Timoshenko beam element [26], Mindlin plate element [27], and Kirchhoff plate element [28]. However, some complex structures [29] in MEMS, such as twisted microplates, curved microbeams, and microsprings, cannot be numerically calculated by these elements used in specific situations, or they can be calculated, but only rough results can be obtained. Therefore, a solid element that exhibits universality may compensate for the shortcomings of some special elements. In this regard, three solid elements, namely, an 8-node element with full integration (FI), an 8-node element with reduced integration (RI), and a 20-node element with FI, are constructed to analyze the scale effect of a generalized elastic body.In the current study, a finite element equation for a generalized elastic body is established on the basis of the principle of virtual work. Furthermore, a penalty function term is introduced into the virtual work equation. The generalized elastic body was made discrete via an 8-node linear element and a 20-node quadratic element. In the numerical examples, the case of a cantilever microbeam was used to discuss the convergence and locking problem in solid elements. The influence of the penalty parameter on the results was also analyzed. Lastly, the scale effects of the deflections of the cantilevered microbeams and the simply supported circular microplates and the frequencies corresponding to the modes of the pretwisted plate were numerically analyzed.


## 2. Generalized Elasticity Theory

A displacement gradient is used to characterize infinitesimal deformation, which can be divided into a symmetric strain tensor and an antisymmetric rotation tensor, i.e.,(1)u∇=ε+Ω.The strain tensorε and the rotation tensor Ω can be written in the following forms:(2)ε=12u∇+∇u,(3)Ω=12u∇−∇u,where u and ∇ are the displacement vector and Laplace operator, respectively.An antisymmetric rotation tensor only has three independent variables, and thus, it can be equivalent to a rotation vectorω, i.e.,(4)ω=−12Ω:∈=12∇×u,where ∈ denotes the permutation tensor.Assuming that there is a line element dx in an elastic body, the increment of the line element is represented by(5)du=ε+Ω⋅dx=ε⋅dx+ω×dx.The tension-compression for a line element represents the translational deformation described by the strain tensor, whereas the bending for the line element represents the rotational deformation measured by the rotation tensor.The degree of rotational deformation is characterized by the curvature tensor(6)χ=∇ω=12∇∇×u.The boundary conditions for a generalized elastic body subjected to body couplec, surface couple m¯, body force b, and surface force t¯ can be expressed as(7)n⋅t=t¯,n⋅m=m¯,onSσ,u=u¯,ω=ω¯,onSu,where n, ω¯, t, and m are the normal vector, rotation vector, Cauchy stress tensor, and couple stress tensor at the boundary surface, respectively.The equations for momentum and moment of momentum of a generalized elastic body are(8)∇⋅t+ρb=ρu¨,∇⋅m+t:∈+ρc=0,where ρ is the density and u¨ is the acceleration vector.In the framework of GET, the Cauchy stress is asymmetrical rather than symmetrical. Therefore, it can be divided into a symmetric stress tensor and an antisymmetric stress tensor, i.e.,(9)t=σ+τ,where the symmetric stress σ follows the generalized Hooke’s law in CET, and the antisymmetric stress τ can be expressed in the form of(10)τ=−12∇⋅m+ρc⋅∈.The constitutive equations for GET are written as(11)σ=2με+λtrεI,(12)m=4ηχ,where λ and μ are Lamé parameters, I is a second-order identity tensor, and trε represents the trace of the strain tensor. η = μl2 is the rotation modulus, where l is the length scale parameter. The constitutive relation between symmetric stress and symmetric strain describes the characteristic of translational deformation, and that between couple stress and curvature defines the feature of rotational deformation. When l = 0, the GET degenerates into the CET.


## 3. Finite Element Equation

The finite element equation for a generalized elastic body is established by applying the principle of virtual work. For a linear elastic body with volumeV and surface area S, the total virtual work done by the external force on the virtual displacement and rotation is written in the following form:(13)δWext=∫Vρδu⋅b+δω⋅cdV+∫Sδu⋅t¯+δω⋅m¯dS−∫Vρδu⋅u¨dV.Meanwhile, the virtual work done by the internal force is written as(14)δWint=∫Vδε:σ+δχ:mdV.By combining equations (13) and (14), the virtual work equation is obtained as(15)∫Vδε:σ+δχ:mdV=∫Vρδu⋅b+δω⋅cdV+∫Sδu⋅t¯+δω⋅m¯dS−∫Vρδu⋅u¨dV.The curvature in equation (15) is the second derivative of the displacement. If the displacement is used as the independent variable, the displacement and its first derivative need to satisfy the continuity between the elements. In order to reduce the continuity requirement, the penalty function method in the constrained variational principle is used here; that is, the displacement and rotation of the nodes are used as independent variables to weaken the continuity between the elements. This means that the continuity requirements for the displacement and its first derivative are converted to continuity requirements for the displacement and rotation. Simultaneously, the rotation φ at the nodes and the rotation vector ω determined by the displacement satisfy the additional condition φ-ω = 0.The additional condition is introduced into the virtual work equation (15) in the form of a product by using a penalty function. Furthermore, by replacing the rotation vector ω in equation (15) with the rotation φ at the nodes, a constrained variational equation can be written as(16)∫Vδε:σ+δχ:mdV−∫Vρδu⋅b+δφ⋅cdV−∫Sδu⋅t¯+δφ⋅m¯dS+∫Vρδu⋅u¨dV+α∫Vδφ−ω⋅φ−ωdV=0,where α is the penalty parameter. The approximate solution obtained by the constrained variational equation approximately satisfies the additional condition, where the larger the α is, the better satisfied the additional condition is.The 8-node hexahedral isoparametric element and the 20-node hexahedral isoparametric element are used to make the generalized elastic body discrete, as shown in Figure1.Figure 1
8-node linear element and 20-node quadratic element.The shape functions [30] of the two elements are as follows:Linear element (8 nodes):(17)Nm=181+ξmξ1+ηmη1+ζmζ,wherem=1,2,3,…,8.Quadratic element (20 nodes):(18)Nn=18ξn2ηn2ζn21+ξnξ1+ηnη1+ζnζξnξ+ηnη+ζnζ−2+14ηn2ζn21−ξ21+ηnη1+ζnζ1−ξn2+14ζn2ξn21−η21+ζnζ1+ξnξ1−ηn2+14ξn2ηn21−ζ21+ξnξ1+ηnη1−ζn2,wheren=1,2,3,…,20.The generalized displacement of the element can be written as(19)u˜=N˜1,N˜2,…,N˜td˜1,d˜2,…,d˜tT=N˜d˜,where N˜k=NkI300NkI3 and d˜k=uφT=ukxukyukzφkxφkyφkzT are the generalized interpolation function and generalized displacement of the kth node in an element, respectively. Moreover, t is the total number of nodes of an element (t = 8 for the linear element and t = 20 for the quadratic element), and I3 is a 3 × 3 identity matrix.The couple stress tensor, curvature tensor, stress tensor, and strain tensor can be written in the matrix form as follows:(20)m=mxxmyymzzmxymyzmxzmyxmzymzxT,χ=χxxχyyχzzχxyχyzχxzχyxχzyχzxT,σ=σxxσyyσzzσxyσyzσxzT,ε=εxxεyyεzzεxyεyzεxzT.Combining equations (2), (3), (11), and (12), the expansion matrix of strain and curvature and that of stress and couple stress can be written as(21)ε˜=εχ=Lu00Lφu˜=Lu˜=LN˜d˜=B˜d˜,(22)σ˜=σm=Du00Dφεχ=D˜ε˜,where B˜ and D˜ are the generalized strain matrix and generalized elasticity matrix, respectively. Du and Dφ=4ηI9×9 are the elasticity matrix in CET and the constitutive relation matrix between couple stress and curvature, respectively. L, Lu, and Lφ are differential operator matrices, and Lu and Lφ are expressed in the form of(23)Lu=∂∂x00∂∂y0∂∂z0∂∂y0∂∂x∂∂z000∂∂z0∂∂y∂∂xT,Lφ=∂∂x00000∂∂y0∂∂z0∂∂y0∂∂x000∂∂z000∂∂z0∂∂y∂∂x000T.Similarly, in consideration of equation (6), φ−ω in equation (16) can be written in the following form:(24)φ−ω=Lαu˜=LαN˜d˜=B˜αd˜,where(25)Lα=0∂2∂z−∂2∂y100−∂2∂z0∂2∂x010∂2∂y−∂2∂x0001is the differential operator matrix of the penalty term and B˜α is the strain matrix of the penalty term.The matrix form of equations (16) can be expressed as(26)∫Vδε˜Tσ˜dV−∫Vρδu˜Tb˜dV−∫Sδu˜Tt˜dS+∫Vρδu˜Tu˜¨dV+α∫Vδφ−ωTφ−ωdV=0.Furthermore, substituting equations (21), (22), and (24) into equation 26 yields(27)δd˜Tρ∫VN˜TN˜dVd˜¨+∫VΒ˜TD˜Β˜dV+α∫VΒ˜αTΒ˜αdVd˜−∫VρN˜Tb˜dV−∫sN˜Tt˜dS=0,where b˜=bcT and t˜=t¯m¯T are the generalized body-force matrix and the generalized surface-force matrix, respectively.Given the arbitrariness of variationδd˜, the finite element formulation can be obtained from equation (27), i.e.,(28)Med˜¨e+Ke+αKαed˜e=Pfe+Pse.For equation (28), the element mass matrix Me, element stiffness matrix Ke, element penalty term matrix Kαe, element body-force matrix Pfe, and element surface-force matrix Pse can be solved using Gaussian integrals, i.e.,(29)Me=ρ∫VN˜TN˜dV=ρ∑k=1r∑j=1r∑i=1rHiHjHkN˜TN˜J,Ke=∫VΒ˜TD˜Β˜dV=∑k=1r∑j=1r∑i=1rHiHjHkΒ˜TD˜Β˜J,Kαe=∫VΒ˜αTΒ˜αdV=∑k=1s∑j=1s∑i=1sHiHjHkΒ˜αTΒ˜αJ,Pfe=∫VρN˜Tb˜dV=ρ∑k=1r∑j=1r∑i=1rHiHjHkN˜Tb˜J,Pse=∫sN˜Tt˜dS=∑k=1r∑j=1rHkHjAN˜Tt˜,t˜acts on the surface ofξ=1,where J is the Jacobian matrix, r and s are the numbers of Gaussian integration points, and Hi, Hj, and Hk are the weight coefficients that correspond to Gaussian integration points. A represents the transformation variable of the load surface between local coordinates and global coordinates, and its three components are(30)Aξ=1=∂y∂η∂z∂ζ−∂y∂ζ∂z∂η2+∂z∂η∂x∂ζ−∂z∂ζ∂x∂η2+∂x∂η∂y∂ζ−∂x∂ζ∂y∂η21/2,Aη=1=∂y∂ζ∂z∂ξ−∂y∂ξ∂z∂ζ2+∂z∂ζ∂x∂ξ−∂z∂ξ∂x∂ζ2+∂x∂ζ∂y∂ξ−∂x∂ξ∂y∂ζ21/2,Aζ=1=∂y∂ξ∂z∂η−∂y∂η∂z∂ξ2+∂z∂ξ∂x∂η−∂z∂η∂x∂ξ2+∂x∂ξ∂y∂η−∂x∂η∂y∂ξ21/2.The structural finite element equation is obtained as(31)Md¨+K+αKαd=P,where M, K, Kα, and P are the assembled mass, stiffness, penalty term stiffness, and load matrices obtained by the usual addition of element matrices.When calculating natural frequency, the characteristic equation is expressed as(32)K+αKα−ω2M=0,where ω represents the natural circular frequency.For the static problem, the mass matrix in equation (31) is omitted and the structural finite element equation is rewritten as(33)K+αKαd=P.


## 4. Numerical Discussions

### 4.1. Influence of Penalty Parameter and Integration Schemes

A cantilever microbeam that is subjected to a uniform downward surface load ofP = 50 MPa at its free end surface is shown in Figure 2. The material and geometric parameters are provided in Table 1. The cantilever microbeam is numerically simulated via an 8-node FI element, an 8-node RI element, and a 20-node FI element. Then, its deflection is evaluated. Given that the exact solution cannot be obtained, a convergent solution is performed. The deflection at the free end center is plotted over the number of DOFs in Figure 3 for the three elements. The results calculated using the 8-node element with RI and the 20-node element with FI converge to a constant value of 2.17 μm, with the latter converging faster than the former. By contrast, the result of the 8-node FI element does not converge; that is, the locking phenomenon occurred.Figure 2
Cantilever microbeam subjected to the uniform surface loadP.Table 1
Numerical simulation parameters.
E
v
ρ
l
B
H
L
Elastic modulus
Poisson’s ratio
Density
Length parameter
Section width
Section height
Length
209 GPa
0.31
7800 kg/m3
6.01μm
10μm
10μm
120μmFigure 3
Maximum deflection versus the number of DOFs.The above was an analysis of the bending deformation of the cantilever microbeam. Next, the tensile deformation of the cantilever microbeam was analyzed. The surface forceP acting on the free end surface of the cantilever microbeam was replaced as the tensile load T = 5000 MPa in the y-direction, and the displacement at the free end of the cantilever microbeam was calculated.The results almost collectively converge to a constant value of 2.84μm, as shown in Figure 4. The convergence of the three elements under bending and tension differs because the linear element with FI cannot withstand bending but can withstand tension. If the linear element with FI is used, then the stiffness will increase and the deflection will decrease (Figure 3). However, the RI scheme may compensate for the problem of excessive stiffness and achieve a better result. For the quadratic element, the results exhibit higher accuracy and faster convergence than the results of the linear element with RI because the former can withstand bending and has a higher-order shape function.Figure 4
Axial displacement versus the number of DOFs.The penalty parameter also affects the calculation results in addition to the DOF. The assessment of the effect of the DOF on the convergence is based on the fact that the effect of the penalty parameter on the results has converged. Here, we check the influence of the penalty parameters. At the same time, to avoid the influence of the DOF on the calculation result, the simulation parameters of the above examples and the maximum DOF when converging are still used here. The deflection at the free end of the cantilever microbeam under bending load is calculated with different penalty parameters. The result is shown in Figure5. Similar to the convergence result of the refined mesh (Figure 3), the 8-node RI and 20-node FI elements exhibit convergence, whereas the 8-node FI element demonstrates divergence.Figure 5
Deflection at the free end of the cantilever microbeam versus the penalty parameter.In addition, the influence of the penalty parameter on the cantilever microbeam under tension is also tested. The axial displacement at the free end of the cantilever microbeam is plotted in Figure6. The axial displacement calculated by the three elements converges when the penalty parameter is less than 1011. However, when the penalty parameter is greater than 1011, the 8-node FI element diverges, whereas the 8-node RI and 20-node FI elements continue to converge, although their values decrease slightly.Figure 6
Axial displacement at the free end of the cantilever microbeam versus the penalty parameter.In general, a principle in which the penalty parameter is sufficiently large should be followed [30]. However, when the penalty parameter approaches infinity, equation (32) degenerates to(34)Kαd≈Pα⟶0.Therefore, the penalty stiffness matrix must remain singular to obtain a nonzero solution. The penalty stiffness matrices calculated using the 8-node RI and 20-node FI elements are singular matrices, whereas that calculated using the 8-node FI element is a nonsingular matrix. Thus, only the 8-node FI element exhibits divergence (Figure5). In addition, the penalty term is constructed via the rotation vector and rotation at the nodes, which indicates that the penalty term is related to the rotational deformation. The cantilever microbeam under tension yields translational deformation that is independent of the rotational deformation, so the penalty term is nearly ineffective for the translational deformation of the cantilever microbeam. However, the degenerated equation (34) appears when the penalty parameter continues to increase, which will lead to divergence of the 8-node FI element because the penalty stiffness matrix calculated using this element is nonsingular. This also explains why the solutions shown in Figure 6 for the three elements converge when the penalty parameter is less than 1011, and the solution of the 8-node FI element diverges when the penalty parameter is greater than 1011.In general, a quadratic element can obtain results with higher accuracy and faster convergence. However, its computational cost is higher than that of a linear element. Therefore, selecting suitable elements is critical for ensuring high efficiency and cost ratio. However, the 8-node FI element should be clearly avoided as much as possible in numerical calculations.


### 4.2. Scale Effect of a Cantilever Microbeam and Simply Supported Microplate

To illustrate the universality of solid elements, a slender cantilever microbeam with a length-to-thickness ratio of 12, a short cantilever microbeam with a length-to-thickness ratio of 6, a thin circular microplate with a diameter-to-thickness ratio of 100, and a medium-thick circular microplate with a diameter-to-thickness ratio of 20 are considered. On the basis of time consumption and calculation accuracy, the 8-node RI element is used to simulate the cantilever microbeam problem, whereas the 20-node FI element is used to analyze the simply supported microplate problem.A uniform downward load ofP = 5000 MPa acts on the free end of the cantilever microbeam and the center of the simply supported circular microplate, as shown in Figures 2 and 7. The geometric parameters of the cantilever microbeam are as follows: width B = 10 μm, height H = 10 μm, and length L = 60 μm and 120 μm. Meanwhile, the geometric parameters of the simply supported circular microplate are as follows: circular plate diameter D = 600 μm, load surface diameter d = 40 μm, and plate thickness Hc = 6 μm and 30 μm. The material parameters of the microbeam and microplate are provided in Table 1. The numerical results calculated in accordance with CET and GET are plotted in Figures 8 and 9, respectively. Notably, Figure 9 adopts a double Y-axis system. That is, the left axis corresponds to the solid and dashed lines, whereas the right axis corresponds to the asterisk and plus signs.Figure 7
Schematic diagram of the simply supported circular plate.Figure 8
Deflection of the cantilever microbeam.Figure 9
Deflection of the simply supported circular microplate.The deflection of the cantilever microbeam in they-direction is shown in Figure 8, whereas the deflection of the simply supported microplate along the radial direction is given in Figure 9. These figures clearly present that the deflections calculated in accordance with GET are smaller than those predicted in accordance with CET.


### 4.3. Scale Effect of Natural Frequencies

The literature on the scale effect of the natural frequency of elastic bodies at the microscale is abundant. However, the object of such studies is mostly regular plates or beams, and only the scale effect of fundamental frequencies is involved. Therefore, the scale effect of frequencies that correspond to different modes of irregular structures should be an interesting research issue. In this section, the scale effect of the natural frequencies that correspond to the modes of the pretwisted plate is numerically analyzed using the 8-node element with RI. The free end of the pretwisted plate (lengthL' = 150 μm, width B' = 30 μm, and thickness h = 3 μm) has a torsional angle ψ = 30° relative to the fixed end, as shown in Figure 10. The material parameters are consistent with those of the cantilever beam in the previous section. The ratio of the thickness of the pretwisted plate to the length scale parameter is regarded as the embodiment of dimensional change, and the geometric dimension ratio of the pretwisted plate remains fixed. The normalized frequencies (i.e., the ratio of the frequency calculated in accordance with GET to the frequency calculated in accordance with CET) that correspond to the three modes are plotted in Figure 11.Figure 10
Configuration of a pretwisted plate. (a) Top view (ψ=0). (b) Main view.
(a)
(b)Figure 11
Curves of normalized frequencies withh/l.As shown in Figure11, the normalized frequency that corresponds to the tension-compression mode remains constant, whereas the frequencies that correspond to the torsion and bending modes gradually increase with a decrease in h/l, and the frequency that corresponds to the torsion mode is higher than the frequency that corresponds to the bending mode. The reason for this phenomenon is that the tension-compression mode involves a translational deformation that is independent of rotational deformation, and therefore, the frequency that corresponds to the tension-compression mode exerts no scale effect. The bending mode and torsion mode involve rotational deformation, which is amplified at the microscale, thereby resulting in the scale effect of the frequencies that correspond to the bending mode and torsion mode. However, no evident scale effect is observed at the macroscale because rotational deformation at the macroscale is minimal and nearly ineffective (Figure 11).In summary, not all frequencies exert a scale effect, and only the frequency corresponding to the mode with rotational deformation has a scale effect.


## 4.1. Influence of Penalty Parameter and Integration Schemes

A cantilever microbeam that is subjected to a uniform downward surface load ofP = 50 MPa at its free end surface is shown in Figure 2. The material and geometric parameters are provided in Table 1. The cantilever microbeam is numerically simulated via an 8-node FI element, an 8-node RI element, and a 20-node FI element. Then, its deflection is evaluated. Given that the exact solution cannot be obtained, a convergent solution is performed. The deflection at the free end center is plotted over the number of DOFs in Figure 3 for the three elements. The results calculated using the 8-node element with RI and the 20-node element with FI converge to a constant value of 2.17 μm, with the latter converging faster than the former. By contrast, the result of the 8-node FI element does not converge; that is, the locking phenomenon occurred.Figure 2
Cantilever microbeam subjected to the uniform surface loadP.Table 1
Numerical simulation parameters.
E
v
ρ
l
B
H
L
Elastic modulus
Poisson’s ratio
Density
Length parameter
Section width
Section height
Length
209 GPa
0.31
7800 kg/m3
6.01μm
10μm
10μm
120μmFigure 3
Maximum deflection versus the number of DOFs.The above was an analysis of the bending deformation of the cantilever microbeam. Next, the tensile deformation of the cantilever microbeam was analyzed. The surface forceP acting on the free end surface of the cantilever microbeam was replaced as the tensile load T = 5000 MPa in the y-direction, and the displacement at the free end of the cantilever microbeam was calculated.The results almost collectively converge to a constant value of 2.84μm, as shown in Figure 4. The convergence of the three elements under bending and tension differs because the linear element with FI cannot withstand bending but can withstand tension. If the linear element with FI is used, then the stiffness will increase and the deflection will decrease (Figure 3). However, the RI scheme may compensate for the problem of excessive stiffness and achieve a better result. For the quadratic element, the results exhibit higher accuracy and faster convergence than the results of the linear element with RI because the former can withstand bending and has a higher-order shape function.Figure 4
Axial displacement versus the number of DOFs.The penalty parameter also affects the calculation results in addition to the DOF. The assessment of the effect of the DOF on the convergence is based on the fact that the effect of the penalty parameter on the results has converged. Here, we check the influence of the penalty parameters. At the same time, to avoid the influence of the DOF on the calculation result, the simulation parameters of the above examples and the maximum DOF when converging are still used here. The deflection at the free end of the cantilever microbeam under bending load is calculated with different penalty parameters. The result is shown in Figure5. Similar to the convergence result of the refined mesh (Figure 3), the 8-node RI and 20-node FI elements exhibit convergence, whereas the 8-node FI element demonstrates divergence.Figure 5
Deflection at the free end of the cantilever microbeam versus the penalty parameter.In addition, the influence of the penalty parameter on the cantilever microbeam under tension is also tested. The axial displacement at the free end of the cantilever microbeam is plotted in Figure6. The axial displacement calculated by the three elements converges when the penalty parameter is less than 1011. However, when the penalty parameter is greater than 1011, the 8-node FI element diverges, whereas the 8-node RI and 20-node FI elements continue to converge, although their values decrease slightly.Figure 6
Axial displacement at the free end of the cantilever microbeam versus the penalty parameter.In general, a principle in which the penalty parameter is sufficiently large should be followed [30]. However, when the penalty parameter approaches infinity, equation (32) degenerates to(34)Kαd≈Pα⟶0.Therefore, the penalty stiffness matrix must remain singular to obtain a nonzero solution. The penalty stiffness matrices calculated using the 8-node RI and 20-node FI elements are singular matrices, whereas that calculated using the 8-node FI element is a nonsingular matrix. Thus, only the 8-node FI element exhibits divergence (Figure5). In addition, the penalty term is constructed via the rotation vector and rotation at the nodes, which indicates that the penalty term is related to the rotational deformation. The cantilever microbeam under tension yields translational deformation that is independent of the rotational deformation, so the penalty term is nearly ineffective for the translational deformation of the cantilever microbeam. However, the degenerated equation (34) appears when the penalty parameter continues to increase, which will lead to divergence of the 8-node FI element because the penalty stiffness matrix calculated using this element is nonsingular. This also explains why the solutions shown in Figure 6 for the three elements converge when the penalty parameter is less than 1011, and the solution of the 8-node FI element diverges when the penalty parameter is greater than 1011.In general, a quadratic element can obtain results with higher accuracy and faster convergence. However, its computational cost is higher than that of a linear element. Therefore, selecting suitable elements is critical for ensuring high efficiency and cost ratio. However, the 8-node FI element should be clearly avoided as much as possible in numerical calculations.


## 4.2. Scale Effect of a Cantilever Microbeam and Simply Supported Microplate

To illustrate the universality of solid elements, a slender cantilever microbeam with a length-to-thickness ratio of 12, a short cantilever microbeam with a length-to-thickness ratio of 6, a thin circular microplate with a diameter-to-thickness ratio of 100, and a medium-thick circular microplate with a diameter-to-thickness ratio of 20 are considered. On the basis of time consumption and calculation accuracy, the 8-node RI element is used to simulate the cantilever microbeam problem, whereas the 20-node FI element is used to analyze the simply supported microplate problem.A uniform downward load ofP = 5000 MPa acts on the free end of the cantilever microbeam and the center of the simply supported circular microplate, as shown in Figures 2 and 7. The geometric parameters of the cantilever microbeam are as follows: width B = 10 μm, height H = 10 μm, and length L = 60 μm and 120 μm. Meanwhile, the geometric parameters of the simply supported circular microplate are as follows: circular plate diameter D = 600 μm, load surface diameter d = 40 μm, and plate thickness Hc = 6 μm and 30 μm. The material parameters of the microbeam and microplate are provided in Table 1. The numerical results calculated in accordance with CET and GET are plotted in Figures 8 and 9, respectively. Notably, Figure 9 adopts a double Y-axis system. That is, the left axis corresponds to the solid and dashed lines, whereas the right axis corresponds to the asterisk and plus signs.Figure 7
Schematic diagram of the simply supported circular plate.Figure 8
Deflection of the cantilever microbeam.Figure 9
Deflection of the simply supported circular microplate.The deflection of the cantilever microbeam in they-direction is shown in Figure 8, whereas the deflection of the simply supported microplate along the radial direction is given in Figure 9. These figures clearly present that the deflections calculated in accordance with GET are smaller than those predicted in accordance with CET.


## 4.3. Scale Effect of Natural Frequencies

The literature on the scale effect of the natural frequency of elastic bodies at the microscale is abundant. However, the object of such studies is mostly regular plates or beams, and only the scale effect of fundamental frequencies is involved. Therefore, the scale effect of frequencies that correspond to different modes of irregular structures should be an interesting research issue. In this section, the scale effect of the natural frequencies that correspond to the modes of the pretwisted plate is numerically analyzed using the 8-node element with RI. The free end of the pretwisted plate (lengthL' = 150 μm, width B' = 30 μm, and thickness h = 3 μm) has a torsional angle ψ = 30° relative to the fixed end, as shown in Figure 10. The material parameters are consistent with those of the cantilever beam in the previous section. The ratio of the thickness of the pretwisted plate to the length scale parameter is regarded as the embodiment of dimensional change, and the geometric dimension ratio of the pretwisted plate remains fixed. The normalized frequencies (i.e., the ratio of the frequency calculated in accordance with GET to the frequency calculated in accordance with CET) that correspond to the three modes are plotted in Figure 11.Figure 10
Configuration of a pretwisted plate. (a) Top view (ψ=0). (b) Main view.
(a)
(b)Figure 11
Curves of normalized frequencies withh/l.As shown in Figure11, the normalized frequency that corresponds to the tension-compression mode remains constant, whereas the frequencies that correspond to the torsion and bending modes gradually increase with a decrease in h/l, and the frequency that corresponds to the torsion mode is higher than the frequency that corresponds to the bending mode. The reason for this phenomenon is that the tension-compression mode involves a translational deformation that is independent of rotational deformation, and therefore, the frequency that corresponds to the tension-compression mode exerts no scale effect. The bending mode and torsion mode involve rotational deformation, which is amplified at the microscale, thereby resulting in the scale effect of the frequencies that correspond to the bending mode and torsion mode. However, no evident scale effect is observed at the macroscale because rotational deformation at the macroscale is minimal and nearly ineffective (Figure 11).In summary, not all frequencies exert a scale effect, and only the frequency corresponding to the mode with rotational deformation has a scale effect.


## 5. Conclusions

The deformation of an elastic body is further classified into translational deformation and rotational deformation. Strain measures the degree of translational deformation, whereas curvature measures the degree of rotational deformation. Together, translational and rotational deformations constitute the generalized elasticity model, which is a supplement to CET. The essence of the scale effect is that the rotational deformation of an elastic body at the microscale is amplified, thereby resulting in different outcomes from CET. However, rotational deformation at the macroscale is extremely small and nearly ineffective, and thus, the scale effect can hardly be observed at the macroscale.The finite element equation of a generalized elastic body is established by following the principle of virtual work, and the penalty function method in the constrained variational principle is adopted by imposing rotation and displacement as independent variables. The convergence of the three elements and the influence of the penalty parameter on the results are discussed. The results show that the 8-node RI and 20-node FI elements exhibit better performance than the 8-node FI element, which demonstrates poor performance and thus should be abandoned.The analysis of the scale effect of the cantilever microbeam and simply supported microplate shows that the 8-node RI and 20-node FI elements can reflect the universality of the solid elements and accurately measure the scale effect. In the analysis of the dynamic characteristics of the pretwisted plate, the frequencies that correspond to the bending and torsion modes exhibit scale effect, whereas the frequency that corresponds to the tension-compression mode does not demonstrate the scale effect. The numerical examples indicate that rotational deformation, rather than translational deformation, causes the scale effect.



---

*Source: 2307689-2019-12-07.xml*