Research Article | Open Access
Numerical Research on a Three-Dimensional Solid Element Based on Generalized Elasticity Theory
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.
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) . Subsequently, Mindlin developed a simpler constitutive model with only one additional material parameter . However, these two constitutive models do not address the issue of the symmetry of the couple stress tensor. Yang et al.  argued that the couple stress tensor is symmetric and proposed a new constitutive model with one additional material parameter. By contrast, Hadjesfandiari and Dargush  believed that the couple stress tensor is antisymmetric and presented a new constitutive model with one additional material parameter. Huang and Huang  and Liu et al.  asserted that the couple stress tensor is asymmetric, and the latter renamed Mindlin’s couple stress theory  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 , Timoshenko beam element , Mindlin plate element , and Kirchhoff plate element . However, some complex structures  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.,
The strain tensor and the rotation tensor can be written in the following forms:where 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.,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
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
The boundary conditions for a generalized elastic body subjected to body couple , surface couple , body force , and surface force can be expressed aswhere n, , , and 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 arewhere ρ is the density and 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.,where the symmetric stress follows the generalized Hooke’s law in CET, and the antisymmetric stress can be expressed in the form of
The constitutive equations for GET are written aswhere λ and μ are Lamé parameters, I is a second-order identity tensor, and 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 volume V 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:
Meanwhile, the virtual work done by the internal force is written as
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 aswhere α 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 Figure 1.
The shape functions  of the two elements are as follows: Linear element (8 nodes): Quadratic element (20 nodes):
The generalized displacement of the element can be written aswhere and 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 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:
Combining equations (2), (3), (11), and (12), the expansion matrix of strain and curvature and that of stress and couple stress can be written aswhere and are the generalized strain matrix and generalized elasticity matrix, respectively. and are the elasticity matrix in CET and the constitutive relation matrix between couple stress and curvature, respectively. , , and are differential operator matrices, and and are expressed in the form of
Similarly, in consideration of equation (6), in equation (16) can be written in the following form:whereis the differential operator matrix of the penalty term and is the strain matrix of the penalty term.
The matrix form of equations (16) can be expressed as
Given the arbitrariness of variation , the finite element formulation can be obtained from equation (27), i.e.,
For equation (28), the element mass matrix , element stiffness matrix , element penalty term matrix , element body-force matrix , and element surface-force matrix can be solved using Gaussian integrals, i.e.,where 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
The structural finite element equation is obtained aswhere 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 aswhere ω 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
4. Numerical Discussions
4.1. Influence of Penalty Parameter and Integration Schemes
A cantilever microbeam that is subjected to a uniform downward surface load of P = 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.
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 force P 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.
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 Figure 5. 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.
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 Figure 6. 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.
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 (Figure 5). 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 of P = 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.
The deflection of the cantilever microbeam in the y-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 (length L' = 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.
As shown in Figure 11, 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.
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.
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by NSAF (Grant no. U1830115) and NSFC (Grant nos. 11772071 and 11372365).
- M. H. Korayem and A. H. Korayem, “Modeling of AFM with a piezoelectric layer based on the modified couple stress theory with geometric discontinuities,” Applied Mathematical Modelling, vol. 45, pp. 439–456, 2017.
- W. Yang, D. He, and W. Chen, “A size-dependent zigzag model for composite laminated micro beams based on a modified couple stress theory,” Composite Structures, vol. 179, pp. 646–654, 2017.
- Y. Sun, J. Cheng, Z. Wang, Y. Yu, L. Tian, and J. Lu, “Analytical approximate solution for nonlinear behavior of cantilever FGM MEMS beam with thermal and size dependency,” Mathematical Problems in Engineering, vol. 2019, Article ID 9637048, 10 pages, 2019.
- M. Ghavami, S. Azizi, and M. R. Ghazavi, “On the dynamics of a capacitive electret-based micro-cantilever for energy harvesting,” Energy, vol. 153, pp. 967–976, 2018.
- S. Jabbari Behrouz, O. Rahmani, and S. Amirhossein Hosseini, “On nonlinear forced vibration of nano cantilever-based biosensor via couple stress theory,” Mechanical Systems and Signal Processing, vol. 128, pp. 19–36, 2019.
- M. H. Ghayesh and H. Farokhi, “Nonsymmetric nonlinear dynamics of piezoelectrically actuated beams,” Journal of Vibration and Acoustics, vol. 141, no. 5, Article ID 051012, 11 pages, 2019.
- H. Salehipour, H. Nahvi, and A. R. Shahidi, “Exact closed-form free vibration analysis for functionally graded micro/nano plates based on modified couple stress and three-dimensional elasticity theories,” Composite Structures, vol. 124, pp. 283–291, 2015.
- A. R. Askari and M. Tahani, “Size-dependent dynamic pull-in analysis of geometric non-linear micro-plates based on the modified couple stress theory,” Physica E: Low-Dimensional Systems and Nanostructures, vol. 86, pp. 262–274, 2017.
- V. S. Khorasani and M. Bayat, “Bending analysis of FG plates using a general third-order plate theory with modified couple stress effect and MLPG method,” Engineering Analysis with Boundary Elements, vol. 94, pp. 159–171, 2018.
- S. Soni, N. K. Jain, and P. V. Joshi, “Analytical modeling for nonlinear vibration analysis of partially cracked thin magneto-electro-elastic plate coupled with fluid,” Nonlinear Dynamics, vol. 90, no. 1, pp. 137–170, 2017.
- S. Soni, N. K. Jain, and P. V. Joshi, “Vibration and deflection analysis of thin cracked and submerged orthotropic plate under thermal environment using strain gradient theory,” Nonlinear Dynamics, vol. 96, no. 2, pp. 1575–1604, 2019.
- S. Soni, N. K. Jain, and P. V. Joshi, “Stability and dynamic analysis of partially cracked thin orthotropic microplates under thermal environment: an analytical approach,” Mechanics Based Design of Structures and Machines, 2019.
- A. W. McFarland and J. S. Colton, “Role of material microstructure in plate stiffness with relevance to microcantilever sensors,” Journal of Micromechanics and Microengineering, vol. 15, no. 5, pp. 1060–1067, 2005.
- S. Hassanpour and G. R. Heppler, “Micropolar elasticity theory: a survey of linear isotropic equations, representative notations, and experimental investigations,” Mathematics and Mechanics of Solids, vol. 22, pp. 224–242, 2015.
- R. D. Mindlin and H. F. Tiersten, “Effects of couple-stresses in linear elasticity,” Archive for Rational Mechanics and Analysis, vol. 11, no. 1, pp. 415–448, 1962.
- R. D. Mindlin, “Influence of couple-stresses on stress concentrations,” Experimental Mechanics, vol. 3, no. 1, pp. 307-308, 1963.
- F. Yang, A. C. M. Chong, D. C. C. Lam, and P. Tong, “Couple stress based strain gradient theory for elasticity,” International Journal of Solids and Structures, vol. 39, no. 10, pp. 2731–2743, 2002.
- A. R. Hadjesfandiari and G. F. Dargush, “Couple stress theory for solids,” International Journal of Solids and Structures, vol. 48, no. 18, pp. 2496–2510, 2011.
- K. Z. Huang and Y. G. Huang, Solid Constitutive Relation, Tsinghua University Press, Beijing, China, 1999, in Chinese.
- Z. Liu, S. Yan, and Z. Fu, “Dynamic analysis on generalized linear elastic body subjected to large scale rigid rotations,” Applied Mathematics and Mechanics, vol. 34, pp. 1001–1016, 2013.
- S. Chakravarty, A. R. Hadjesfandiari, and G. F. Dargush, “A penalty-based finite element framework for couple stress elasticity,” Finite Elements in Analysis and Design, vol. 130, pp. 65–79, 2017.
- B. T. Darrall, G. F. Dargush, and A. R. Hadjesfandiari, “Finite element Lagrange multiplier formulation for size-dependent skew-symmetric couple-stress planar elasticity,” Acta Mechanica, vol. 225, no. 1, pp. 195–212, 2014.
- X. Ma and W. Chen, “Refined 18-DOF triangular hybrid stress element for couple stress theory,” Finite Elements in Analysis and Design, vol. 75, pp. 8–18, 2013.
- X. Ma and W. Chen, “24-DOF quadrilateral hybrid stress element for couple stress theory,” Computational Mechanics, vol. 53, no. 1, pp. 159–172, 2014.
- M. Kandaz and H. Dal, “A comparative study of modified strain gradient theory and modified couple stress theory for gold microbeams,” Archive of Applied Mechanics, vol. 88, no. 11, pp. 2051–2070, 2018.
- A. T. Karttunen, J. Romanoff, and J. N. Reddy, “Exact microstructure-dependent Timoshenko beam element,” International Journal of Mechanical Sciences, vol. 111-112, no. 112, pp. 35–42, 2016.
- B. Zhang, Y. He, D. Liu, Z. Gan, and L. Shen, “A non-classical Mindlin plate finite element based on a modified couple stress theory,” European Journal of Mechanics—A/Solids, vol. 42, pp. 63–80, 2013.
- M. Tahani, A. R. Askari, Y. Mohandes, and B. Hassani, “Size-dependent free vibration analysis of electrostatically pre-deformed rectangular micro-plates based on the modified couple stress theory,” International Journal of Mechanical Sciences, vol. 94-95, no. 95, pp. 185–198, 2015.
- N. Lobontiu and E. Garcia, Mechanics of Microelectromechanical Systems, Kluwer Academic Publishers, New York, NY, USA, 2005.
- O. Zienkiewicz, R. Taylor, and J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, Beijing Word Publishing Corporation, Beijing, China, 7th edition, 2015.
Copyright © 2019 Zhen Qiao and Zhanfang Liu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.