A New Nonlinear Spatial Compliance Model Method for Flexure Leaf Springs with Large Width-to-Length Ratio under Large Deformation

Flexure leaf spring (FLS) with large deformation is the basic unit of compliant mechanisms with large stroke. The stiffness along the non-working directions of FLSs with large width-to-length ratio (w/L) is high. The motion stability of the compliant mechanism based on this type of FLS is high. When this type of FLS is loaded along the width direction, the shear deformation needs to be characterized. Nevertheless, currently available compliance modeling methods for FLS are established based on Euler–Bernoulli beam model and cannot be used to characterize shear models. Therefore, these methods are not applicable in this case. In this paper, a new six-DOF compliance model for FLSs with large w/L is established under large deformation. The shear deformation along the width direction model is characterized based on the Timoshenko beam theory. The new constraint model and differential equations are established to obtain a high-precision compliance model expression for this type of FLS. The effects of structural parameters on the compliance of the FLS are analyzed. Finally, the accuracy of the model is verified both experimentally and by finite element simulation. The relative error between theoretical result and experiment result is less than 5%.


Introduction
Compliant mechanisms are mechanical mechanisms that rely on the elastic deformation of materials and, among other uses, have been widely applied in medical devices [1][2][3], precision instruments [4][5][6], sensors [7][8][9], and microelectronic-mechanical systems (MEMS) [10,11]. Compliant mechanisms can be divided into centralized and distributed compliant mechanisms. The deformation of centralized compliant mechanisms is only concentrated at flexure hinges under small deformation. The compliance of the flexure hinges and the overall centralized compliant mechanism can be modelled based on the premise of linear deformation [12][13][14]. On the other hand, distributed compliant mechanisms are based on the deformation of the entire structure or most of it. In the same design space, the stroke of the distributed compliant mechanism is larger than that of the centralized one; thus, distributed compliant mechanisms have been widely used in precision machinery and precision instruments [15][16][17]. An accurate compliance model is the basis for the structural design of such mechanisms. When large deformation occurs in a distributed compliant mechanism, there is geometric nonlinearity in the flexure leaf spring (FLS). To this end, it is important to investigate and develop a nonlinear compliance model with high accuracy able to simulate the large deformation of FLS. transformation method. In Section 4, the effects of the structural parameters on FLS compliance and compliance ratio are analyzed. In Section 5, an experimental platform is built to verify the accuracy of the model proposed in this paper.

Expressions of Deformation and Load
The deformation of the FLS under load and the definition of the structural parameters are given in Figure 1. The position of the point P1 on the FLS is expressed as (X, 0, 0) when the FLS is not deformed and as P 1 (X + U X , U Y , U Z ) after deformation. It is defined by the cross-section where the point P 1 (point P 1 after deformation) is located perpendicular to the tangential direction of the FLS axis at P 1 after deformation. A coordinate system is established on the cross-section, where the point P 1 is located after deformation. The X d1 axis of the coordinate system coincides with the tangent of the axis at P 1 , and when the FLS is not deformed, the directions of the Y d1 and Z d1 axes are the same as those of the Y and Z axes of the initial coordinate system. The transformation relationship between the initial coordinate system (X, Y, Z) and the one on the cross-section where the point P 1 is located after deformation (X d1 , Y d1 , Z d1 ) is depicted in Figure 2 and it can be expressed as: where c(Θ) = cos(Θ), s(Θ) = sin(Θ).
Micromachines 2022, 13, x FOR PEER REVIEW 4 of 27    When the deformation of the FLS is small, i.e., less than 0.1 L, the following approximate relationship can be satisfied [31]: The torsion around the X axis and the curvature around the Y and Z axes can be expressed as [28]: The loads on the end of the FLS are F XL , F YL , F ZL , M XL , M YL , and M ZL . The deformation and rotation angles at the point PL at the end of the FLS are expressed as U X1 , U Y1 , U Z1 , and Θ Xd1 , Θ Y1 , Θ Z1 , respectively. The load applied to any point on the FLS can be expressed as (4).

Differential Governing Equations
This paper considers that the w/L of the FLS is large. This section analyzes the calculation process of the deformation of every DOF.

Deformation in the Z-O-X Plane
The loads in the Z-O-X plane are F ZL and M YL . Due to the F ZL , there is shear deformation in the FLS width direction, while, due to the M YL , there is deflection deformation. The deformation in the Z-O-X plane can be expressed based on the Timoshenko beam theory as follows: where k 2 is the shear coefficient, here k 2 = 0.85; E is the elastic modulus; G is the shear modulus; A is cross-sectional area of the FLS, which is defined as A = wt; and I Y is the inertia moment around the Y axis, which is defined as I Y = w 3 t/12.

Deformation in the Y-O-X Plane
The loads in the Y-O-X plane are F YL and M ZL , due to which there is deflection deformation in the FLS thickness direction. The deformation in the Y-O-X plane can be expressed based through the Euler-Bernoulli beam theory as follows: where I Z is the inertia moment around the Z axis, which is defined as I Z = wt 3 /12.

Torsion and Extension/Contraction
When the FLS is loaded by M X , there is warping deformation due to the effect of shear force [30]. The deformation can be expressed by the following equation: where J is the torsional constant of the FLS and can be expressed as [34]: (8) where η = t/w. The expression of f(η) is as follows: The relationship between the expansion/contraction deformation of the FLS and F X can be expressed as follows [31]: It should be noted that κ 2 X is so small that it can be ignored.

Global System of Equations
The global system of differential governing equations for the spatial FLS is as follows: Micromachines 2022, 13, 1090 6 of 25 The load and deformation can be normalized as follows: The system of equations can be expressed as:

Differential Governing Equations
In Figure 1, the FLS is fixed at O; thus, the boundary conditions for the displacement of the FLS can be expressed as: Every deformation of the FLS is derivable from the axial coordinate X and can be expressed as a power series of X. Therefore, the solution method based on series is employed to solve the system of equations. Every deformation can be expressed as: where x = X/L. The expressions of θ Y and θ Z can be obtained according to Equation (2) as follows: The following system of equations can be obtained by substituting Equation (16) into Equation (11): The coefficients in Equation (16) are related to the load and structural parameters of the FLS. When Equation (18) is solved, the coefficients of the same order terms in the system of equations are made be equal to zero. In this paper, the numerical solution method is applied. The a n , b n , c n , and d n coefficients gradually approach 0 with increasing n. Therefore, the result can be considered to be convergent. When the deformation of the FLS is small, the high-order terms in Equation (16) are smaller than the first-and second-order terms. Consequently, the fifth-and higher-order terms are ignored. The final result is obtained using the polynomial fitting method and is expressed as Equation (19).

Chain Model of the Spatial FLS
The chain model of a spatial FLS with large deformation is illustrated in Figure 3. The forces and torques loaded on the end of the FLS are F xO , F yO , F zO and M xO , M yO , M zO , respectively, while the displacement and rotation angles of the FLS end are p N , q N , r N and θ xdN , θ yN , θ zN , respectively. Subsequently, the FLS is divided into several flexure elements with equal length L i . Every flexure element has two endpoints and is modeled with the six-DOF compliance model introduced in Section 2. The load and deformation are defined in Figure 4. It is assumed that the FLS is divided into N flexure elements. The local coordinate system x i -y i -z i is established at the endpoint (node i) of the element i and along its tangent direction. The local coordinate system of the first element is established at the fixed end of the FLS (node 0). sponding displacements and rotation angles are Δpi, Δqi, Δri and Δθxdi, Δθyi, Δθzi, respectively. The transformation matrix between the coordinate system on the end of the i-th flexure element and that on the endpoint (node i) is given as Equation (27):     The transformation matrix between the local coordinate system xi-yi-zi of the i-th element and the global coordinate system X-Y-Z is as follows: which can be defined as Equation (29): sponding displacements and rotation angles are Δpi, Δqi, Δri and Δθxdi, Δθyi, Δθzi, respectively. The transformation matrix between the coordinate system on the end of the i-th flexure element and that on the endpoint (node i) is given as Equation (27):      The transformation matrix between the local coordinate system xi-yi-zi of the i-th element and the global coordinate system X-Y-Z is as follows: which can be defined as Equation (29): The free end of the FLS is defined as node N. The forces and torques loaded on the i-th flexure element are ∆F xi , ∆F yi , ∆F zi and ∆M xi , ∆M yi , ∆M zi , respectively, while the corresponding displacements and rotation angles are ∆p i , ∆q i , ∆r i and ∆θ xdi , ∆θ yi , ∆θ zi , respectively. The transformation matrix between the coordinate system on the end of the i-th flexure element and that on the endpoint (node i) is given as Equation (27): The transformation matrix between the local coordinate system x i -y i -z i of the i-th element and the global coordinate system X-Y-Z is as follows: which can be defined as Equation (29): The value of every angle can be calculated using the above transformation matrix, and the calculation process is given as Equation (30): Equation (32) can be deduced: By further transforming Equations (32) and (33) can be derived as follows: The total displacement of node i in the chain model can be expressed as follows: The loads applied on node i are: The transformation relationship between the loads in the local coordinate system x i -y i -z i and those in the global coordinate system X-Y-Z is: Furthermore, the relationship between the forces exerted on adjacent nodes is: The chain model of the spatial FLS can be obtained by combining the equation systems of every flexure element and using the six-DOF compliance model introduced in Section 2 to model every flexure element. Subsequently, a relationship between loads (F xO , F yO , and F zO ), torques (M xO , M yO , and M zO ), displacements (p N , q N , and r N ), and rotation angles (θxdN, θyN, and θzN) can be obtained. Due to the complexity of the equation system, Newton's method is used to solve it. The calculation results (p N , q N , r N and θ xdN , θ yN , θ zN ) gradually approach a fixed value with an increasing number of flexure elements n. Consequently, the chain model can be considered as convergent.

Compliance
According to the previous sections, the relationship between deformation and load of the spatial FLS under large deformation is nonlinear. In this section, the effects of the structural parameters on the compliance of the FLS as a function of the deformation along the thickness direction (working direction) are analyzed. Based on the variables in the chain model introduced in Section 2.3, the compliance in every direction can be defined as follows: The range of each parameter is as follows: 1.
The range of t is [0.6 mm, 1 mm].
The material of the FLS was set as spring steel (60Si2Mn) with a Young's modulus of E = 206 GPa and Poisson's ratio of v = 0.29. The effects of the different parameters are presented in Figures 5-8, where (a-c) show the effects of L, w, and t, respectively, on compliance.   In Figure 5, it can be observed that the value of cy decreased with the increase of the FLS deformation along the thickness direction. Moreover, the value of cy increased with the increase of L, decreased with the increase of w, and increased with the increase of t. The variation range of cy increased with increasing L, decreased with increasing w, and decreased with increasing t over the deformation range of the FLS.
As can be seen in Figure 6, the value of cz increased with the increase of the FLS deformation along the thickness direction. The value of cz increased with increasing L, decreased with increasing w, and increased with increasing t. The variation range of cz increased with increasing L, decreased with increasing w, and decreased with increasing t over the deformation range of the FLS.
According to Figure 7, the value of cx increased with the increase of the FLS deformation along the thickness direction. The value of cx increased with increasing L, decreased with increasing w, and decreased with increasing t. The variation range of cx increased with increasing L, decreased with increasing w, and decreased with increasing t over the deformation range of the FLS.
In Figure 8, it can be seen that the value of crx decreased with the increase of the FLS deformation along the thickness direction. More specifically, the value of crx increased with increasing L, decreased with increasing w, and decreased with increasing t. The variation range of crx did not change with the variation of L and that of t but decreased with increasing w over the deformation range of the FLS. In Figure 5, it can be observed that the value of c y decreased with the increase of the FLS deformation along the thickness direction. Moreover, the value of c y increased with the increase of L, decreased with the increase of w, and increased with the increase of t. The variation range of c y increased with increasing L, decreased with increasing w, and decreased with increasing t over the deformation range of the FLS.
As can be seen in Figure 6, the value of c z increased with the increase of the FLS deformation along the thickness direction. The value of c z increased with increasing L, decreased with increasing w, and increased with increasing t. The variation range of c z increased with increasing L, decreased with increasing w, and decreased with increasing t over the deformation range of the FLS.
According to Figure 7, the value of c x increased with the increase of the FLS deformation along the thickness direction. The value of c x increased with increasing L, decreased with increasing w, and decreased with increasing t. The variation range of c x increased with increasing L, decreased with increasing w, and decreased with increasing t over the deformation range of the FLS.
In Figure 8, it can be seen that the value of c rx decreased with the increase of the FLS deformation along the thickness direction. More specifically, the value of c rx increased with increasing L, decreased with increasing w, and decreased with increasing t. The variation range of c rx did not change with the variation of L and that of t but decreased with increasing w over the deformation range of the FLS.

Compliance Ratio
The ratio of the compliance along the working direction to that along the non-working directions is an important indicator that can reflect the motion stability of the FLS. The compliance ratio indicates whether the performance of the FLS can ensure the smoothness of motion in the working direction and resist the disturbance forces in the non-working directions. To this end, the following compliance ratios were analyzed: The range of each parameter is as follows: 1.
The range of t is [0.6 mm, 1 mm].
The material of the FLS was set as spring steel (60Si2Mn) with a Young's modulus of E = 206 GPa and Poisson's ratio of v = 0.29. The effects of the different parameters are exhibited in Figures 9-11, where (a-c) show the effects of L, w, and t, respectively, on compliance ratio.

Compliance Ratio
The ratio of the compliance along the working direction to that along the non-working directions is an important indicator that can reflect the motion stability of the FLS. The compliance ratio indicates whether the performance of the FLS can ensure the smoothness of motion in the working direction and resist the disturbance forces in the non-working directions. To this end, the following compliance ratios were analyzed: The range of each parameter is as follows:   As can be seen in Figure 9, the value of c z /c y increased with the increase of the FLS deformation along the thickness direction. The value of c z /c y did not change with the variation of L, decreased with increasing w and increased with increasing t. The variation range of c z /c y did not change with the variation of L, decreased with increasing w, and increased with increasing t over the deformation range of the FLS.
In Figure 10, it can be observed that the value of c x /c y increased with the increase of the FLS deformation along the thickness direction. The value of c x /c y decreased with increasing L, remained unchanged with the variation of w, and increased with increasing t. The variation range of c x /c y decreased with increasing L, did not change with the variation of w, and increased with increasing t over the deformation range of the FLS.
According to Figure 11, the value of c rx /c y decreased first and then increased with the increase of the FLS deformation along the thickness direction. The value of c rx /c y decreased with increasing L, remained unchanged with the variation of w, and increased with increasing t. The variation range of c rx /c y decreased with increasing L, did not change with the variation of w, and increased with increasing t over the deformation range of the FLS.

Theoretical Model Verification
In this paper, the compliance model of the spatial large deformation FLS is verified both experimentally and through finite element simulations. The structure of the FLS, the nodes of the chain model, and the positions of the applied load are depicted in Figure 12. The structural parameters of the FLS are given in Table 2. The chain model described in Section 3 was used to discretize the FLS. More specifically, the FLS was divided into 18 flexure elements; the length of each element was 10 mm. The holes for applying different loads to the FLS were placed on the boundary where node 16 was located. The distance between adjacent holes was 20 mm. In both the simulation and experiment, the deformation and rotation angle of the midpoints (P 1 , P 2 , P 3 ) on the discrete boundary where nodes 3, 8, and 13 are located were compared with the calculation results obtained by the theoretical model. In order to analyze the rotation angle of the FLS in the simulation and experimental results, the deformations of the points (P 11 , P 12 , P 21 , P 22 , P 31 , P 32 ) on the boundary line on the FLS where nodes 3, 8, and 13 are located were extracted. The distance between the different points is given in Table 3. The coordinates of the points when the FLS was not deformed are listed in Table 4. The coordinates of each point on the FLS after deformation are listed in Table 5 and are expressed based on the chain model introduced in Section 3. experimental results, the deformations of the points (P11, P12, P21, P22, P31, P32) on the bound ary line on the FLS where nodes 3, 8, and 13 are located were extracted. The distance be tween the different points is given in Table 3. The coordinates of the points when the FLS was not deformed are listed in Table 4. The coordinates of each point on the FLS after deformation are listed in Table 5 and are expressed based on the chain model introduced in Section 3.  Nodes in the chain model (P 1 , P 2 , P 3 ) and nodes used for displacement analysis (P 11 , P 12 , P 21 , P 22 , P 31 , P 32 ). (0.03, 0, −0.03) P 12 (0.03, 0, 0.03) P 2 (0.08, 0, 0.) P 21 (0.08, 0, −0.03) P 22 (0.08, 0, 0.03) P 3 (0.13, 0, 0.) P 31 (0.13, 0, −0.03) P 32 (0.13, 0, 0.03) Table 5. Coordinates of points after deformation.

Finite Element Simulation
The mesh of the FLS model is depicted in Figure 13. The model was meshed with tetrahedral elements, the number of nodes was 228,711 and the number of elements was 111,623. In this paper, the deformation of the FLS under loads applied at three different positions was simulated. The applied load was 1 N, and the application positions of the loads are given in Table 6. The simulation results are demonstrated in Figure 14, and the points P 11 , P 12 , P 21 , P 22 , P 31 , P 32 were recorded from the simulation results. The displacement and rotation angle of each point (P 1 , P 2 , P 3 ) could be calculated based on the simulation results of points (P 11 , P 12 , P 21 , P 22 , P 31 , P 32 ). The theoretical calculation and simulation results of points (P 1 , P 2 , P 3 ) are listed in Table 7. It can be found that the error was smaller than 5%. This verifies that the theoretical model proposed in this paper is accurate.

Experiment
The experimental platform used to test the FLS is depicted in Figure 15. A 6-DOF manipulator (UNIVERSAL ROBOTS, UR5) carrying a surface structured light sensor (TECHLEGO, Q3; measurement accuracy of ±0.005 mm) was used to scan the deformed FLS and capture its spatial contours. The spatial change of the position of points (P11, P12, P21, P22, P31, P32) on the FLS after deformation was obtained by processing the spatial contour data. Subsequently, the deformation and rotation angle of points (P1, P2, P3) on the FLS were calculated. The FLS used in the experiment is depicted in Figure 16

Experiment
The experimental platform used to test the FLS is depicted in Figure 15. A 6-DOF manipulator (UNIVERSAL ROBOTS, UR5) carrying a surface structured light sensor (TECHLEGO, Q3; measurement accuracy of ±0.005 mm) was used to scan the deformed FLS and capture its spatial contours. The spatial change of the position of points (P11, P12, P21, P22, P31, P32) on the FLS after deformation was obtained by processing the spatial contour data. Subsequently, the deformation and rotation angle of points (P1, P2, P3) on the FLS were calculated. The FLS used in the experiment is depicted in Figure 16

Experiment
The experimental platform used to test the FLS is depicted in Figure 15. A 6-DOF manipulator (UNIVERSAL ROBOTS, UR5) carrying a surface structured light sensor (TECHLEGO, Q3; measurement accuracy of ±0.005 mm) was used to scan the deformed FLS and capture its spatial contours. The spatial change of the position of points (P11, P12, P21, P22, P31, P32) on the FLS after deformation was obtained by processing the spatial contour data. Subsequently, the deformation and rotation angle of points (P1, P2, P3) on the FLS were calculated. The FLS used in the experiment is depicted in Figure 16. The FLS

Experiment
The experimental platform used to test the FLS is depicted in Figure 15. A 6-DOF manipulator (UNIVERSAL ROBOTS, UR5) carrying a surface structured light sensor (TECH-LEGO, Q3; measurement accuracy of ±0.005 mm) was used to scan the deformed FLS and capture its spatial contours. The spatial change of the position of points (P 11 , P 12 , P 21 , P 22 , P 31 , P 32 ) on the FLS after deformation was obtained by processing the spatial contour data. Subsequently, the deformation and rotation angle of points (P 1 , P 2 , P 3 ) on the FLS were calculated. The FLS used in the experiment is depicted in Figure 16. The FLS was marked at points (P 11 , P 12 , P 21 , P 22 , P 31 , P 32 ) in order to facilitate the subsequent data processing. Due to that the scanning area of the structured light sensor was small and could not fully cover the FLS, the FLS needed to be marked at reference points to implement the splicing algorithm of spatial contour. Then, the scanned contour data of the entire FLS could be obtained. The FLS deformed under load is displayed in Figure 17. The load was m = 100 g and the gravitational acceleration g = 9.8066 N/kg. The spatial contours of the deformed FLS are presented in Figure 18. The experimental results of the points (P 1 , P 2 , P 3 ) were obtained after processing the spatial contour data and are listed in Table 8. The error was smaller than 5%. Consequently, the accuracy of the theoretical model proposed in this paper was again verified. was marked at points (P11, P12, P21, P22, P31, P32) in order to facilitate the subsequent data processing. Due to that the scanning area of the structured light sensor was small and could not fully cover the FLS, the FLS needed to be marked at reference points to implement the splicing algorithm of spatial contour. Then, the scanned contour data of the entire FLS could be obtained. The FLS deformed under load is displayed in Figure 17. The load was m = 100 g and the gravitational acceleration g = 9.8066 N/kg. The spatial contours of the deformed FLS are presented in Figure 18. The experimental results of the points (P1, P2, P3) were obtained after processing the spatial contour data and are listed in Table 8.
The error was smaller than 5%. Consequently, the accuracy of the theoretical model proposed in this paper was again verified.   was marked at points (P11, P12, P21, P22, P31, P32) in order to facilitate the subsequent data processing. Due to that the scanning area of the structured light sensor was small and could not fully cover the FLS, the FLS needed to be marked at reference points to implement the splicing algorithm of spatial contour. Then, the scanned contour data of the entire FLS could be obtained. The FLS deformed under load is displayed in Figure 17. The load was m = 100 g and the gravitational acceleration g = 9.8066 N/kg. The spatial contours of the deformed FLS are presented in Figure 18. The experimental results of the points (P1, P2, P3) were obtained after processing the spatial contour data and are listed in Table 8.
The error was smaller than 5%. Consequently, the accuracy of the theoretical model proposed in this paper was again verified.

Conclusions
In this paper, a 6-DOF compliance model for FLSs with large w/L under large deformation has been proposed. A new differential equation system is established based on Timoshenko beam theory along the width direction. The shear deformation along the width direction of the FLS with large w/L and the spatial geometric nonlinearity can be accurately described by this method. In Section 2, a compliance model for flexure elements with small deformation (i.e., <0.1 L) has been proposed. Then, the chain model of the spatial FLS was established in Section 3. The flexure elements with small deformation are assembled and expanded to accurately represent the large deformation of the FLS. Based on the chain model and the compliance model of the flexure elements, the effects of the structural parameters on the compliance and the compliance ratio of the FLS under large deformation were analyzed in Section 4. In Section 5, the theoretical model was verified experimentally and through finite element simulations, and the accuracy of the proposed theoretical model was proven. It was found that the relative error between theoretical result and experiment result was less than 5%.