Abstract
The guide stiffness performance directly affects the motion of the micromechanism in accuracy and security. Therefore, it is crucial to analyze the guide stiffness precisely. In this paper, a high-precision guide stiffness analysis method for the micromechanism by the boundary element method (BEM) is proposed. The validity and accuracy of the analysis method are tested by a guide stiffness experiment. In order to ensure the accuracy and safety during the micromechanism motion, a guiding unit of the micromechanism was designed based on the guiding principle. The guiding unit can provide parasitic motion and additional force in the motion of the micromechanism. Then, the stiffness equations of the beam element are derived by the boundary element method. The stiffness equation of straight circular flexure hinge is analyzed by rigid discretization and rigid combination, and the guide stiffness of the mechanism is investigated by rigid combination. Finally, according to the actual situation, the stiffness matrix of the guide rail (Kb) was proposed, and the analytical value of the guide stiffness was calculated to be 22.2 N/μm. The guide stiffness performance experiment was completed, and the experimental value is 22.3 N/μm. Therefore, the error between the analysis method and the experimental results is 0.45%. This study provides a new method for the stiffness analysis of high-precision micromechanisms and presents a reference for the design and stiffness analysis of complex structures. This method is helpful for stiffness analysis of the microrotary mechanism with high accuracy.
1. Introduction
With the rapid development of numerical control technology, advanced manufacturing has called for higher requirements for precision machining technology [1–3]. Macro/micro-dual-drive technology has solved the contradiction of large-motion travel and high-motion precision for a single-drive mechanism system and provided motion for large travel and high precision. The macro/micro-dual-drive system has been widely used in high-technology fields, such as the advanced weapons industry, biological medicine, and precision electronics [4–9]. In the macro/micro-dual-drive system, the microdrive system can compensate the motion error of the macrodrive system and provide high-precision motion for the macro/micro-dual-drive system [10–13]. As this mechanism can provide accurate motion based on the compound motion principle of a flexure hinge, the micromechanism is a core component of the microdrive system [14–17], and the performance of the micromechanism is directly affected by the performance of the microdrive system [18–21].
Many scholars have conducted research on the motion performance and practical application of the microlinear mechanism [22–25]. Meanwhile, some scholars focus on studying microrotary mechanism [26, 27]. Because of parasitic motion and additional force, the motion accuracy and safety of the micromechanism will be reduced. The parasitic movement in nonmotion directions will affect the accuracy of the movement in motion direction. Microactuators (i.e., PZT) are generally assembled inside the micromechanism to drive the mechanism, and usually, the stiffness of microactuators in motion direction is very high. The stiffness of microactuators in nonmotion directions is usually very low, and additional force will badly affect the safety of microactuators during motion of a micromechanism. As the guiding motion principle, the guide unit of micromechanism can provide relaxation of parasitic movement and additional force during the motion. Therefore, it is necessary to design a guide unit of micromechanism.
As a significant performance indicator of the micromechanism, the stiffness can affect the motion performance of the mechanism under external load. In recent years, many scholars have conducted research on the stiffness of the micromechanism, and significant achievements have been achieved. Using the relational expression of Castigliano’s second theorem to express the displacement and force of the end of the flexure hinge, the relative displacement vector of the flexure hinge was analyzed and calculated, and the stiffness matrix of the system was proposed, and meanwhile, the minimum motion principle was obtained [28]. The stiffness model is based on the way the flexure members are connected in serial or parallel combinations, and the modeling allows one to formulate the functional relationship between stiffness and dimensions as well as the free shape of the FPM in the design process, while the proposed analytical model is validated by the FEM model and experiments [29]. A novel six-strut compliant parallel mechanism based on the development of wide-range flexure hinges is explored, the stiffness equation of an individual flexure hinge is established, and the stiffness of the whole mechanism is modeled by assembling stiffness matrixes and formulating constraint equations, after which the system stiffness influence plots are presented and discussed [30]. A bridge micromechanism with 3 degrees of freedom was designed. The z-direction of the mechanism had higher rigidity and stability, which could resist large external thrust. The structure of the organization is optimized by adopting the sequential quadratic rule. Finally, the experiments show that the optimization method is highly effective [31]. A novel type of flexible mechanism composed of two straight circular flexure hinges with multiple notches was designed. The topology method is used to optimize the mechanism structure, and the finite element method was used to analyze the stiffness, rotation accuracy, and horizontal stress of the mechanism. The analysis results show that the mechanism performances are excellent [32]. A fully compliant, potentially monolithic, and power transmission mechanism that can rectify a large lateral offset between two parallel rotational axes is presented. The transmission stiffness and the actuation stiffness of the designed device are predicted by the theoretical model and finite element modeling [33].
As the guide unit can ensure accuracy and safety during the motion of the micromechanism, it is significant to analyze the guide stiffness of the micromechanism. However, there are rarely relevant papers to explore the guide stiffness of the micromechanism. The stiffness of the guide unit is interpreted by the stiffness formula in most studies. This method will gain an inevitable error because of the stiffness model error in the analysis of a micromechanism with a complex structure. BEM is a method of dividing discrete elements within a defined range on the boundary, accurately decomposing the research object and then analyzing the boundary problem. BEM has many advantages, such as less preparation work in data, being vital to deal with particular engineering problems, and high accuracy [34, 35]. Moreover, BEM has a significant advantage in solving practical engineering problems with complex structures or complex boundary conditions. It is a particularly prominent advantage of BEM to analyze beams with exact solutions. Therefore, it is appropriate to examine the guide stiffness by BEM.
In this paper, a guide unit of micromechanism is designed for relaxation of parasitic movement and additional force in the micromechanism movement, which can ensure accuracy and safety during its motion, and this research is beneficial to the structural design of the micromechanism. A high-precision guide stiffness analysis method by the boundary element method (BEM) is proposed, and the validity and accuracy of the analysis method are tested by experiment. And, this study provides a new method to analyze the stiffness of micromechanism with high precision.
The rest of this paper is organized as follows. Section 2 introduces the design and guiding motion principle of the guiding unit, which can provide relaxation of parasitic movement and additional force. The high-precision guide stiffness analysis method by the boundary element method (BEM) is given in Section 3. In Section 4, the guide stiffness matrix (Kb) and analytical value are calculated. Section 5 carries out the guide stiffness experiment and analyses results. The conclusions are drawn in Section 6.
2. The Guide Unit of the Micromechanism
When flexible material is stressed, the flexure hinge will produce minor deformation, realizing the transmission motion and guiding motion. In all kinds of flexure hinges, the straight circular flexure hinge is the most widely used because it can meet demand both precision and range of motion.
The equivalent model of straight circular flexure hinge is shown in Figure 1. Point a is the starting point of the hinge, point c is the midpoint of the hinge, and point b is the terminal point of the hinge. Point a′ and point b′ are the starting point and the terminal point the hinge after motioning. li is half the length of the hinge in the x-direction.

(a)

(b)

(c)
The guiding motion principle is a micromechanism motion principle based on symmetrical structure, to ensure micromechanism motion is accurate in the motion direction. The displacement that can provide motion only exists in the motion direction and there is no nonmotion direction during the motion process.
2.1. Design of the Guide Unit
We have designed two kinds of micromechanisms in our research project, as shown in Figure 2. The microlinear mechanism is shown in Figure 2(a), which can realize transmission in the linear motion with different ratios, include amplification and reduction. And, the microrotary mechanism is shown in Figure 2(b), which can transform an input of linear motion into an output of rotary motion according to a particular relationship. For relaxation of parasitic movement and additional force, the guide unit is designed in two kinds of micromechanisms. In Figure 2(a), the mechanism has eight units: 1-2, 3-4, 5-6, 7-8, 9-10, 11-12, 13-14, and 15-16, and they are guide units. Similarly, in Figure 2(b), the mechanism has eight guide units: 6-7, 8-9, 10-11, 12-13, 14-15, 16-17, 18-19, and 20-21.

(a)

(b)
The motion model of the guide unit is shown in Figure 3. The guide unit consists of flexure hinge I, rod 1, and flexure hinge II. The eight guide units are symmetrically distributed on both sides of b parts. When the b parts move along the y-positive direction or the y-positive direction driven by PZT, the eight guide units can be seen as y-direction motion guide rails of b parts by obtaining the same deformations along with two y-directions. Therefore, b parts will bear the load and produce displacement in y-directions and will not bear or make displacement in nonmotion directions (i.e., x-directions). The above motion process can realize the guiding motion of the micromechanism, ensure the accuracy of the motion with no displacement in nonmotion directions, and ensure the safety of PZT working with no load in nonmotion directions.

2.2. Guiding Motion Principle of the Guide Unit
The guide unit is symmetrically distributed on both sides of PZT, which can realize the guiding effect of the micromechanism, and the guide unit is the same as the linear motion guide. The guiding motion principle of the guide unit is shown in Figure 4, 1 is a fixed station, 2 is moving parts, and four guide units (3, 4, 5, 6) constitute the linear motion guide to guiding the linear motion of 2. Figures 4(a) and 4(b), respectively, show the structure diagram of the balance and motion of the guide unit group under the action of driving load P.

(a)

(b)
Taking guide unit 5 as an example, the guiding motion principle of guide unit 5 is shown in Figure 5. The rectangular coordinate system is established, as shown in Figure 5, and the left endpoint of the guide unit 5 is the origin of coordinates. The structural diagram of the guide unit in balance is shown in Figure 5(a), and the right side of the guide unit is fixed. The guide unit consists of flexure hinge I, rod 1 and flexure hinge II. Flexure hinge I is equivalent simplified including revolute pair c and two rods (ac and cb), flexure hinge II is equivalent simplified into revolute pair e and two rods (de and ef); the length of the four rods (ac, cb, de and ef) are li, and the length of rod 1 is l1.

(a)

(b)

(c)
The structural diagram of the guide unit in motion is shown in Figure 5(b). In the left side of the guide unit, a linear motion displacement as under an external load will occur. As the right side of the guide unit is fixed, the guide unit will produce tensile deformation and bending deformation, which is equivalent as all the rods are extended and all the points are rotated. The left side of the guide unit motion is a displacement of in the direction of the y-axis, which can be regarded as the rotation of the guide unit with the direction of . The equivalent figure of guiding motion is shown in Figure 5(c). The guide unit will rotate an angle after the left side of the guide unit moving with under an external load.
The guiding motion of the other three guide units (3, 4, and 6) is similar to that of guiding unit 5. Since the four guide units are symmetrically distributed around the fixed working station 1 and the moving parts 2, the guide unit realizes the motion direction of the guiding moving part 2, and each of the guide units (3, 4, 5, and 6) will rotate an angle after parts 2 moving with .
l is the length of guide unit after motion:
is the included angle between guide unit and x-axis after motion:
is the displacement in y-directions of every position point at guide unit after motion:where is coordinate value in x-directions of every position point at guide unit.
And, equation (3) is the mathematical model of guiding motion.
During the motion of moving parts 2, the guide units will produce force and displacement because of elongation of rods, the forces of the left two guide units (3, 4) can balance the forces of the right two guide units (5, 6). The displacement of the left two guide units (3, 4) can balance the displacements of the right two guide units (5, 6). As the symmetrical structure, force and displacement between the left two guide units (3, 4) and the right two guide units (5, 6) are counteract. This principle can ensure force and displacement of motion only exist in the motion direction and that there is no nonmotion direction force and displacement during motion.
According to the guiding motion principle, the micromechanism has no parasitic motion and additional force in the process of motion, which ensures the accuracy and safety of the movement.
3. The Stiffness Analysis Theory Based on the Boundary Element Method
The guide unit is shown in Figure 3, and it consists of 2 straight circular flexure hinges (i.e., flexure hinge I and flexure hinge II) and a rod (i.e., rod 1), and the structure of the straight circular flexure hinge is more complex than the rod. BEM has advantages such as less preparation work in data, strong capacity to deal with particular engineering problems, and high accuracy, and BEM has a significant advantage in solving practical engineering problems with complex structure and complex boundary conditions. Thus, it is appropriate to analyze the stiffness of the guide unit by BEM.
3.1. Analysis of Static Bending Including Tension and Compression for Beam Element
The static tension and compression for beam elements in the x-y plane are shown in Figure 6. Points 1 and 2 are the two endpoints of the beam, P is the external load of tension and compression for the beam, L is the length of the beam, and the main parameters of the static tension and compression for the beam are shown in Table 1.

The tension and compression stress condition for the microunit of the beam element is shown in Figure 7.

The boundary equation of tension and compression for the beam element can be deduced aswhere is the differential of longitudinal displacement, is the differential of axial force, U is the longitudinal displacement value of the boundary point, and F is the axial force value of the boundary point.
Static bending for the beam element in the x-y plane is shown in Figure 8. Point 1 and point 2 are the two endpoints of the beam, P is the external load of bending for the beam, L is the length of the beam, and the main parameters of the static bending for the beam are shown in Table 2.

The bending stress condition for the microunit of the beam element is shown in Figure 9.

The boundary equation of bending for the beam element can be deduced aswhere q is distributed load, Q is shear force, M is bending moment, E is the longitudinal modulus of elasticity, I is a moment of inertia, V is bending displacement (i.e., deflection), T is bending angular displacement (i.e., deflection angle), Pi is the ith external load, and xi is the position of action for Pi.
The boundary stiffness equation of bending, including tension and compression for the beam element, can be deduced from equations (4) and (5) as
And, equation (6) can be transformed aswhere , , is the bending including tension and compression stiffness of the beam, is the bending including tension and compression displacement of the beam, is the boundary force, and is the external load.
Equation (7) is the boundary stiffness equation of bending, including tension and compression for the beam element.
In equation (7), can be transformed aswhere
And, can be transformed aswherewhere U1 and F1 are the displacements and boundary forces of endpoint 1 and U2 and F2 are the displacement and boundary forces of endpoint 2.
The stiffness matrix [K] can be described as
Therefore, the bending boundary stiffness equation of beam element with tension and compression can be transformed as
The boundary stiffness equation (13) will be used in the rigid combination of the beam element.
3.2. Rigid Combination of the Beam Element
The method of rigid combination is shown in Figure 10. Point 1 and point 2 are the two boundaries of beam a and point 3 and point 4 are the two boundaries of beam b. The beam a and beam b are combined in boundary 3 and boundary 4, as shown in Figure 10(a). The beam c is the new beam through a rigid combination between beam a and beam b, as shown in Figure 10(b), and the beam c can be a rigid combination with other beams.

The stiffness equation of beam a is
And, the stiffness equation of beam b iswhere and are stiffness matrixes of beam a and beam b, and are displacement matrixes of beam a and beam b, and are boundary force matrixes of beam a and beam b, and and are external load matrixes of beam a and beam b.
The force and displacement condition equation of the rigid combination between beam a and beam b is
The rigid combination equation between beam a and beam b with intermediate points iswhere
The rigid combination equation between beam a and beam b without intermediate points iswhere
The stiffness model of the straight circular flexure hinge is a nonuniform beam with a continuously varying height, and the nonuniform beam consists of numerous cross-sections. It can be analyzed by stiffness equations (7) and (13), and rigid combination equations (17) and (19). Then, the guide stiffness matrix can be calculated by a rigid combination between the two straight circular flexure hinges and rod, and the guide stiffness analysis based on the boundary element method is accomplished.
4. Analysis of Guide Stiffness
4.1. Stiffness Calculation of the Straight Circular Flexure Hinge
The stiffness model of the straight circular flexure hinge is a nonuniform beam with a continuously varying height, and the nonuniform beam consists of numerous cross sections. The elastic modulus E is the same between the sections, but section area A and the moment of inertia I are different. If xi is the value of the x-axis for an uncertain ith point i, the variables A and I all have certain function relationships with xi; if α (xi) is the section stiffness function at xi, EI0 is the datum section stiffness, and the section stiffness at xi can be described as
In equation (14), Ii can be solved only in a few cases (such as ). Otherwise, Ii cannot be solved. That is, Ii cannot be solved in any cases by a conventional method. When analyzing the stiffness in any case, the boundary element method is used to discretize the nonuniformity accurately, as shown in Figure 11. Figure 11(a) is the discretization model of a straight circular flexure hinge, Figure 11(b) is the height hi and length Δl of the discretization beam element, and Figure 11(c) is the method of the height hi. The nonuniform beam is discretized of n nonuniform beam elements with equal length. The nonuniform beam element approaches a uniform beam element when the value of n is large enough. We can improve the accuracy of the analysis by increasing the value of n, and hi is the average height of the nonuniform beam element.

(a)

(b)

(c)
The nonuniform beam can be discretized of n uniform beams with equal length and different heights. The height of the uniform beams varies with the xi location variation at the straight circular flexure hinge. The n uniform beams are numbered as shown in Figure 11(a).
If xi is the x-axis value of the endpoint of the ith uniform beam, then xh is the x-axis value of the midpoint of the ith uniform beam.
In order to improve the analysis accuracy, hi is the height of the midpoint location of the ith uniform beam, and hi is
And, the ith uniform beam moment of inertia Ii iswhere R is the radius of the straight circular flexure hinge, t is the minimum height of the straight circular flexure hinge, Δl is the length of each uniform beam (Δl = 2R/n), and xh is the x-axis value of the midpoint of ith uniform beam (xh = xi − 0.5Δl).
Putting the values of Δl (i.e., L) and Ii (i.e., I) and the material and shape parameters of the flexure hinge into equation (7), the stiffness equation of the uniform beam is established. Then, through the rigid combination of uniform beams in equation (13) and computer calculation, the stiffness equation and stiffness matrix of the straight circular flexure hinge are analyzed.
Therefore, the stiffness equation and stiffness matrix Kh of the straight circular flexure hinge can be analyzed as described above.
In Figure 3, the relevant parameters are described as follows: R is 5 mm, t is 2 mm, L (the length of rod 1) is 10 mm, h (the height of rod 1) is 12 mm, b (the thickness of the guide unit) is 12 mm, and the material of the mechanism is 60si2Mn with a longitudinal modulus of elasticity E of 2.06 × 1011 Pa, and n = 1000, while the relevant parameters are set the same as the parameters of the microguide mechanism. Then, the computed result of the stiffness matrix Kh is
In the above calculation process and the matrix Kh, the units are displacements, angular displacements, force, the moment of force in meters (m), degrees (°), Newton (N), and Newton-meters (N·m). The units of relevant parameters of all stiffness matrixes in this paper are defined in the same manner.
4.2. Stiffness Calculation of the Guide Unit
The guide stiffness model of the micromechanism is shown in Figure 3, and the guide unit consists of flexure hinge I, rod 1, and flexure hinge II. The stiffness analysis process of the guide rail is as follows: the rigid combined flexible hinge I and rod 1 are the new rod a, and the rigid combined rod a and the flexible hinge II are the new rod b. Therefore, rod b is the equivalent guide stiffness model of the guide unit. The stiffness matrix Kh of flexure hinge I and flexure hinge II are equal, and the calculation method of Kh is shown in Section 4.1 of this paper.
The stiffness matrix K1 of rod 1 can be calculated by putting the relevant parameters into the matrix K of equation (7) as
The stiffness matrixes of flexure hinge I and rod 1 are Kh and K1, and the stiffness matrix Ka of rod a can be calculated by equation (19) as
The stiffness matrixes of flexure rod a and hinge II are Ka and Kh, and the stiffness matrix Kb of rod b can be calculated by equation (19) as
Hence, the equivalent guide stiffness matrix of the mechanism is stiffness matrix Kb, and matrix Kb is the stiffness matrix of guide stiffness of the micromechanism.
If the left side of the guide unit is fixed and the right side of the guide unit is freed, then the driving force Py = 200 N with a y-positive direction is loaded at the right side of the guide unit, as shown in Figure 3. When guide stiffness of the micromechanism is analyzed if we use the rigid combination equation with intermediate points as equation (17), the bending displacement of deflectionV and the bending angular displacement of deflection angle T at each intermediate point can be calculated. Sixty detection points are set and uniformly distributed in the length direction of the guide unit, and the 60 points are numbered as 1, 2, …, 60 from left to right along the length direction of the guide unit. The deflection curve and the deflection angle curve of the guide unit are shown in Figure 12(a) (the unit of deflections is μm) and Figure 12(b) (the unit of deflection angles is ”).

(a)

(b)
The analysis shows that the right endpoint deflection Δu is 9.0 μm when the driving force Py as 200 N is loaded at the right side of the guide unit, and the stiffness value of the guide unit is K = 22.2 N/μm.
Therefore, the bending stiffness value of the guide unit K is 22.2 N/μm.
5. Experimental Verification
5.1. Experimental Scheme
A guide stiffness experiment has been carried out to verify the effectiveness of the proposed method of guide stiffness analysis. To avoid the influence of external factors such as vibration and temperature change on the experimental results, the experiment was carried out on the vibration isolation platform of the constant temperature laboratory. The guide stiffness experiment is shown in Figure 13, the principle of the guide stiffness experiment is shown in Figure 13(a), and the experiment picture is shown in Figure 13(b). In the guide stiffness experiment, the central experiment test apparatus are microguide mechanism, two displacement sensors are DGS-6C, data acquisition, force sensor is LCM500, force-loading device, and machine tool base.

(a)

(b)
In the experiment, the microguide mechanism is fixed, the different driving force Py is loaded on the right side of the mechanism by a force-loading device, and the driving force Py is detected by a force sensor. The right-side displacement (Δy1) and the left-side displacement (Δy2) of the microguide mechanism are detected by displacement sensor 1 and displacement sensor 2. The right-side displacement (Δy1) is the whole movement displacement of the mechanism driven by an external force. The guide stiffness (K′) of the microguide mechanism in the experiment iswhere ∆u is the right-side relative displacement (relative to the location of the left side) of the microguide mechanism in the y-direction.
5.2. Experimental Results and Analysis
The experimental results are shown in Figure 14. Curve a represents the guide stiffness performance in the loading stage, and the curve of describes the guide stiffness performance in the unloading stage.

For the guide stiffness performance in the loading stage, the linear equation fitted between driving force Py and relative displacement ∆u is
The linearity of the linear equation is 0.9998.
For the guide stiffness performance in the unloading stage, the linear equation fitted between driving force Py and relative displacement ∆u is
The linearity of the linear equation is 0.9997.
As shown in the analysis in Section 4, the guide stiffness (K) calculated by BEM is 22.2 N/μm, and the relationship between K and the two linear equations is shown in Figure 15.

The experimental results show that the guide stiffness ( and ) of the mechanism in the loading stage and unloading stage are 22.1 N/μm and 22.4 N/μm, and the lowest linearity of the bending stiffness linear equation are 0.9998 and 0.9997, which indicates that the experiment results are accurate and believable.
In the loading stage, the experimental guide stiffness = 22.1 N/μm and the calculated guide stiffness K = 22.2 N/μm. Thus, the relative error of the guide stiffness between the two methods is −0.45%. In the unloading stage, the experimental guide stiffness = 22.4 N/μm and the calculated guide stiffness K = 22.2 N/μm. Thus, the relative error of the guide stiffness between the two methods is 0.90%.
Therefore, in the whole stage, the averages experimental guide stiffness = 22.3 N/μm. Thus, the error of the guide stiffness between analysis and experiment is 0.45%. In addition, the experimental results show that the high-precision guide stiffness analysis method of the micromechanism based on BEM is adequate and accurate.
6. Conclusion
This paper designs a micromechanism guiding unit to ensure the accuracy and safety of the micromechanism movement. On this basis, a method for analyzing the stiffness of the micromechanism guide rail based on BEM is proposed. In addition, the guide stiffness experiment was also carried out to test the validity and accuracy of this guide stiffness analysis method.
The guiding unit of the micromechanism is designed for the parasitic movement and the relaxation of the additional force during the movement of the micromechanism to ensure the accuracy and safety of the movement of the micromechanism. This research is conducive to the structural design of the micromechanism. The guide stiffness analysis method of the micromechanism guide rail based on BEM can better solve the special engineering problems. The guide beam stiffness analysis method of the BEM is a high-precision analysis method, and the value of the guide beam stiffness is analytical. The experimental results show that the calculated guide rail stiffness is 22.2 N/μm, the experimental guide rail stiffness is 22.3 N/μm, and the error between the methods is 0.45%. By continuing to study the method (adding rigid combination in the coordinate transformation method), the transformation stiffness of the microstructure can be analyzed. Then, the completion of two kinds of micromechanism (in Section 2.1, as shown in Figure 2) can be undertaken.
This method is suitable for other flexible units and can be used for stiffness analysis of various micromechanisms. Further work will be carried out on the stiffness analysis of the microrotary mechanism with high accuracy by BEM.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no potential conflicts of interest with respect to the research in this article.
Acknowledgments
The authors would like to thank the National Natural Science Foundation of China (nos. 51805428 and 51705417), the Basic Research Project of Natural Science in Shaanxi Province (nos. 2018JQ5205 and 2017JM5042), the Key Research and Development Project of Shaanxi Province (no. 2016KTZDGY4-01), the Science and Technology Plan Project of Xi’an City (no. 201805036YD14CG20(7)), the Special Scientific Research Project of the Shaanxi Province Department of Education (no. 18JK0528), and the Science and Technology Plan Project of Beilin District (no. GX2032) for their support.