Modeling of Motion Characteristics and Performance Analysis of an Ultra-Precision Piezoelectric Inchworm Motor

Ultra-precision piezoelectric inchworm motor (PIM) is widely used in the optical equipment, microelectronics semiconductor industry and precision manufacturing for motion and positioning, but the multi-physics field simulation model for estimating PIM performance and assisting motor design is rarely studied. The simulation model in this paper aimed to provide researchers with direct and convenient PIM performance evaluation to assist the motor design and development. According to the existing advanced inchworm motor products, a multi-physics field coupling model involving solid mechanics and electrostatics using the finite element method (FEM) was established. The motion gesture and performance (driving force and travel) of the PIM were analyzed, respectively. The simulation results showed that the motion gesture of the inchworm motor was well consistent with that of the actual motor product. The driving force from the simulation was close to that of the actual product, and the maximum error was 2.8%. As for the PIM travel, there was a maximum travel error of 0.6 μm between the simulation and official data. The performance parameters of the piezoelectric materials under certain specifications can be simulated by the multi-physics field coupling model. Therefore, the multi-physics field coupling simulation model is suitable for PIM performance evaluation and assisting motor development.


Introduction
A piezoelectric motor is common in ultra-precision motion, positioning, and micro-assembly with its characteristics of high resolution, instant response, and long travel [1][2][3][4]. It supports application scenarios such as the fine tuning of the mirror array of radio telescopes [5], micro adjustments in ultra-precision optical equipment [6], and gesture control of engineering robots [7]. For a piezoelectric motor with an inchworm's walking mode, it is called a PIM. The high-load and ultra-compact PIM has great potential in the field of demanding semiconductor lithography, because it has high requirements for reliability, position resolution, and long-term stability [8]. In particular, the PIM is very likely to be an advanced solution of a chunk motion actuator, replacing the conventional voice coil linear motor in the current lithography chunk.
The working principle of a piezoelectric motor is an inverse piezoelectric effect. That is, when the voltage is applied to an unconstrained piezoelectric material (PM), the PM will produce geometric deformation [9]. The deformation movement is determined by the geometry of the PM, the direction of the electric field, and polarization. For example, the deformation modes of three common piezoelectric ceramics are shown in Figure 1. After polarization, only three piezoelectric strain modes-d 33 , d 15 , and d 31 exist. For the d 33 mode, the deformation direction is consistent with that of the electric field and polarization as shown in Figure 1a; the cuboid structure can also generate the d 33 strain mode. For the d 15 mode, if the electric field direction is on the x-axis of the global coordinate system and the polarization direction is perpendicular to it, the volume deformation (shear deformation) is the vector sum of the PM's deformation on the x-axis and the z-axis as presented in Figure 1b; the deformation exists only on the x-axis and z-axis but not on the y-axis. For the d 31 mode, illustrated in Figure 1c, if the electric field direction is on the z-axis, the deformation direction is on the x-axis. The motion gesture of the inchworm motor in this paper is composed of the d 33 and d 15 modes [10]. The d 33 strain mode is responsible for the extension and contraction of the actuation legs in the vertical direction, while the d 15 strain mode is responsible for generating shear deformation in the horizontal direction, which makes the motor shaft move in the same direction. This direct contact method without a flexible mechanism avoids the mechanical coupling process, so there is no reduction in accuracy and reliability [11].
Materials 2020, 13, x FOR PEER REVIEW  2 of 19 strain modes-d33, d15, and d31-exist. For the d33 mode, the deformation direction is consistent with that of the electric field and polarization as shown in Figure 1a; the cuboid structure can also generate the d33 strain mode. For the d15 mode, if the electric field direction is on the x-axis of the global coordinate system and the polarization direction is perpendicular to it, the volume deformation (shear deformation) is the vector sum of the PM's deformation on the x-axis and the z-axis as presented in Figure 1b; the deformation exists only on the x-axis and z-axis but not on the y-axis. For the d31 mode, illustrated in Figure 1c, if the electric field direction is on the z-axis, the deformation direction is on the x-axis. The motion gesture of the inchworm motor in this paper was composed of the d33 and d15 modes [10]. The d33 strain mode is responsible for the extension and contraction of the actuation legs in the vertical direction, while the d15 strain mode is responsible for generating shear deformation in the horizontal direction, which makes the motor shaft move in the same direction. This direct contact method without a flexible mechanism avoids the mechanical coupling process, so there is no reduction in accuracy and reliability [11].
(a) (b) (c) The key technique of the ultra-precision piezoelectric motor has been widely studied. Simulation modeling is a powerful tool for PIM design in many studies. Li et al. [11] simulated the flexure clamp mechanism + actuator-type PIM using ANSYS FEM to perform stress analysis, but the mechanical coupling process of the PM deformation leading to the movement of the flexure clamp was not involved in the frequency simulation analysis. An ANSYS FEM simulation of a flexure clamping mechanism + middle-drive mechanism-type PIM was conducted by Ma et al. [12], and the related simulation also did not reflect the energy coupling process between the PM and clamping mechanism. Therefore, the common strategy between Li et al. [11] and Ma et al. [12] was to separately simulate the compliant mechanism and the driving mechanism (PM) of the system to indirectly obtain the desired output. As for PIM research in recent years, a four-foot rotary piezoelectric motor without a clamping mechanism was developed in Liu's study [13], its transient simulation work takes into account friction between the rotor and driving feet, but the key deformation mode of the PM and the abovementioned coupling problem was not reflected in the simulation. For the flexure mechanism + drive mechanism-type PIM research, such as Ma et al.'s [12], Zhang et al. [14] conducted stress and frequency analyses of the driving mechanism presented in a steady state simulation. To sum up the simulations in the above studies, the strain mode of the PMs were not involved, the motion gesture of the PIM must be achieved by the correct strain mode of the PM. In addition, the mechanical coupling between the PM and flexure mechanism or motor shaft was rarely reflected. Finally, the travel performance of the PIM was obtained by actual measurement or indirect formula calculation; it was not capable of being obtained directly from the simulation model [12,13,14]. However, for the motor design, the relationship between the specifications of the PM (i.e., type and size) and precision motor performance (i.e., driving force and travel) is what professional developers want to know most. The above cases also appear in References [15,16,17,18,19]. In view of the issues mentioned, a multi-physics field coupling The key technique of the ultra-precision piezoelectric motor has been widely studied. Simulation modeling is a powerful tool for PIM design in many studies. Li et al. [11] simulated the flexure clamp mechanism + actuator-type PIM using ANSYS FEM to perform stress analysis, but the mechanical coupling process of the PM deformation leading to the movement of the flexure clamp was not involved in the frequency simulation analysis. An ANSYS FEM simulation of a flexure clamping mechanism + middle-drive mechanism-type PIM was conducted by Ma et al. [12], and the related simulation also did not reflect the energy coupling process between the PM and clamping mechanism. Therefore, the common strategy between Li et al. [11] and Ma et al. [12] was to separately simulate the compliant mechanism and the driving mechanism (PM) of the system to indirectly obtain the desired output. As for PIM research in recent years, a four-foot rotary piezoelectric motor without a clamping mechanism was developed in Liu's study [13], its transient simulation work took into account friction between the rotor and driving foot, but the key deformation mode of the PM and the abovementioned coupling problem were not reflected in the simulation. For the flexure mechanism + drive mechanism-type PIM research, such as Ma et al.'s [12], Zhang et al. [14] conducted stress and frequency analyses of the driving mechanism presented in a steady state simulation. To sum up the simulations in the above studies, the strain modes of the PM were not involved, the motion gesture of the PIM must be achieved by the correct strain mode of the PM. In addition, the mechanical coupling between the PM and flexure mechanism or motor shaft was rarely reflected. Finally, the travel performance of the PIM was obtained by actual measurement or indirect formula calculation; it was not capable of being obtained directly from the simulation model [12][13][14]. However, for the motor design, the relationship between the specifications of the PM (i.e., type and size) and precision motor performance (i.e., driving force and travel) is what professional developers want to know most. The above cases also appear in References [15][16][17][18][19]. In view of the issues mentioned, a multi-physics field coupling simulation involving electrostatics and solid mechanics using FEM was established for direct PIM performance evaluation to assist the motor design and development under certain PM.
In this paper, according to the advanced inchworm motor, which is developed by PI Motion, the multi-physics coupling simulation was built for PIM performance estimation. During movement, contact friction will occur between the actuation foot and the motor shaft. The contact friction is a highly nonlinear process which undoubtedly increases the complexity of the model and makes the model very difficult to solve. Therefore, the two studies were carried out as follows: (1) focusing on the analysis of the gesture motion of the four-foot PIM, a coupling model without a drive shaft was established; (2) a model with a driving shaft was established to study the motor's performance.

Material
The core piezoelectric material of the PIM is PIC255 which is classification in accordance with EN 50324-1. The PIC255 belongs to "soft" ferroelectric ceramic based on lead zirconate titanate (PZT) and barium titanate. The tetragonal unit cell of "soft" piezo ceramic is shown in Figure 2. The "soft" ideal characteristic of general piezo ceramics can be explained by Figure 3. Figure 3 illustrates the strain-electric field characteristic of "soft" piezo ceramic when a bipolar voltage is applied. When the negative electric field increases to a critical positive electric field, E C , the strain gradually decreases. The strain increases with the increasing electric field when the electric field strength exceeds the E C . If the positive electric field decreases to critical negative electric field −E C , the strain increases gradually. When the electric field is less than −E C , the strain increases with the decreasing electric field. So PIC255 is applied to the actuator's application in dynamic operations (alternating voltage). In addition, the Curie temperature of the PIC255 is 350 • C; it ensures that the PIC255 has extraordinary stability at normal temperature [20]. simulation involving electrostatics and solid mechanics using FEM was established for direct PIM performance evaluation to assist the motor design and development under certain PM. In this paper, according to the advanced inchworm motor, which is developed by PI Motion, the multi-physics coupling simulation was built for PIM performance estimation. During movement, contact friction will occur between the actuation foot and the motor shaft. The contact friction is a highly nonlinear process which undoubtedly increases the complexity of the model and makes the model very difficult to solve. Therefore, the two studies were carried out as follows: (1) focusing on the analysis of the gesture motion of the four-foot PIM, a coupling model without a drive shaft was established; (2) a model with a driving shaft was established to study the motor's performance.

Material
The core piezoelectric material of the PIM is PIC255 which is classification in accordance with EN 50324-1. The PIC255 belongs to "soft" ferroelectric ceramic based on lead zirconate titanate (PZT) and barium titanate. The tetragonal unit cell of "soft" piezo ceramic is shown in Figure 2. The "soft" ideal characteristic of general piezo ceramics can be explained by Figure 3. Figure 3 illustrates the strain-electric field characteristic of "soft" piezo ceramic when a bipolar voltage is applied. When the negative electric field increases to a critical positive electric field, EC, the strain gradually decreases. The strain increases with the increasing electric field when the electric field strength exceeds the EC. If the positive electric field decreases to critical positive electric field −EC, the strain increases gradually. When the electric field is less than −EC, the strain increases with the decreasing electric field. So PIC255 is applied to the actuator's application in dynamic operations (alternating voltage). In addition, the Curie temperature of the PIC255 is 350 °C; it ensures that the PIC255 has extraordinary stability at normal temperature [20].

Simulation
According to PICA Shear actuation and PiezoWalker Motors of PI Motion in References [21,22], the PIM model was built as shown in Figure 4. The system consisted of a fixed axle (domain 1), four actuation legs, and a driving shaft (domain 3). Each actuation leg was composed of upper and lower piezoelectric materials (domain 2.1 and 2.2). In addition, it was necessary to remove the part of the driving shaft in the model when only paying attention to the gesture analysis of the four actuation legs. The three actuator products, P-123.01, P-143.01, and P-153.01, were simulated in this paper [23]. The simulation model was built using COMSOL FEM 5.5 software.

Simulation
According to PICA Shear actuation and PiezoWalker Motors of PI Motion in References [21,22], the PIM model was built as shown in Figure 4. The system consisted of a fixed axle (domain 1), four actuation legs, and a driving shaft (domain 3). Each actuation leg was composed of upper and lower piezoelectric materials (domain 2.1 and 2.2). In addition, it was necessary to remove the part of the driving shaft in the model when only paying attention to the gesture analysis of the four actuation legs. The three actuator products, P-123.01, P-143.01, and P-153.01, were simulated in this paper [23]. The simulation model was built using COMSOL FEM 5.5 software.

Physical Field Analysis
The whole physical model was the coupling of solid mechanics and electrostatics; the model with the driving shaft also involved the highly nonlinear contact friction between the four actuation legs and the driving shaft. Nonlinearity will sharply increase the difficulty of model solving, so the model needs to reduce the number of contact surfaces as much as possible. The minimum number of contact surfaces of the model with a drive shaft is four as presented in Figure 4. This model reduced the nonlinear contact surface as much as possible on the premise of ensuring the integrity of the motion process. As for the model without the driving shaft, the model solution can be easily realized due to the lack of contact friction.

Solid Mechanics
The whole model including domains 1-3 are described by solid mechanics. The two types of piezoelectric constitutive equations in solid mechanics were used to analyze the inverse piezoelectric effect which reflects the coupling between the structural (strain and stress) and electrical domains (the electric field). The strain and stress of the material can be obtained from piezoelectric constitutive equation. The strain and stress of the volume elements on the contact surfaces produce normal and tangential force which is needed in the Coulomb friction model. Finally, the piezoelectric constitutive equation must be applied to all PM regions (domains 2.1 and 2.2).
The strain-charge constitutive equation reflects the coupling relationship between material stress and relative permittivity at constant stress which is shown in Equation (1): where S is the 6 × 1 strain matrix, reflecting the deformation degree of PM. SE is the sixth order flexibility coefficient square matrix, T is the 6 × 1 stress matrix, d is the 3 × 6 piezoelectric strain constant matrix, E is the 3 × 1 electric field strength matrix in the X, Y, and Z direction of the orthogonal global Cartesian coordinate, D is the electric displacement, ε0 is the vacuum dielectric constant, εrT the is relative dielectric constant at constant stress, and t is the transposition symbol. The above coefficient matrix forms are determined by the PIC255 material characteristics. The PIC255 belongs to a tetragonal crystal system, so the flexibility coefficient matrix SE is a symmetrical matrix; it only contains six independent elements as shown in Equation (3):

Physical Field Analysis
The whole physical model was the coupling of solid mechanics and electrostatics; the model with the driving shaft also involved the highly nonlinear contact friction between the four actuation legs and the driving shaft. Nonlinearity will sharply increase the difficulty of model solving, so the model needs to reduce the number of contact surfaces as much as possible. The minimum number of contact surfaces of the model with a drive shaft is four as presented in Figure 4. This model reduced the nonlinear contact surface as much as possible on the premise of ensuring the integrity of the motion process. As for the model without the driving shaft, the model solution can be easily realized due to the lack of contact friction.

Solid Mechanics
The whole model including domain 1, domain 2.1 and domain 2.2, and domain 3 is described by solid mechanics. The two types of piezoelectric constitutive equations in solid mechanics were used to analyze the inverse piezoelectric effect which reflected the coupling between the structural (strain and stress) and electrical domains (the electric field). The strain and stress of the material can be obtained from the piezoelectric constitutive equations. The strain and stress of the volume elements on the contact surfaces produce normal and tangential force which is needed in the Coulomb friction model. Finally, the piezoelectric constitutive equation must be applied to all PM regions (domains 2.1 and 2.2).
The strain-charge constitutive equation reflects the coupling relationship between material stress and relative permittivity at constant stress which is shown in Equation (1): where S is the 6 × 1 strain matrix, reflecting the deformation degree of PM. S E is the sixth order flexibility coefficient square matrix, T is the 6 × 1 stress matrix, d is the 3 × 6 piezoelectric strain constant matrix, E is the 3 × 1 electric field strength matrix in the X, Y, and Z direction of the orthogonal global Cartesian coordinate, D is the electric displacement, ε 0 is the vacuum dielectric constant, ε rT is the relative dielectric constant at constant stress, and t is the transposition symbol. The above coefficient matrix forms are determined by the PIC255 material characteristics. The PIC255 belongs to a tetragonal crystal system, so the flexibility coefficient matrix S E is the symmetrical matrix; it only contains six independent elements as shown in Equation (2): Generally, only three piezoelectric strain constants of the d 33 , d 15 , and d 31 are left after polarization of PM, and its piezoelectric strain coefficient matrix becomes: After PM polarization, the relative permittivity mainly exists in the directions that are both parallel (ε rT33 ) and perpendicular (ε rT11 ) to the polarization direction. The matrix form of the relative permittivity ε rT is as follows: The stress-charge constitutive equation reflects the coupling relationship between material strain and relative permittivity at constant strain. Equation (5) is as below: where c E is the sixth order stiffness coefficient square matrix, e is the 3 × 6 piezoelectric stress constant matrix, also called coupling matrix, and ε rS is the relative dielectric constant at constant strain. The form of the stiffness coefficient matrix c E is the same as that of the flexibility coefficient matrix S E , which is also determined by the PIC 255 crystal system type as presented in Equation (6): Similarly, the form of piezoelectric stress coefficient matrix e strictly follows the form of strain coefficient matrix d; it is listed below: This ε rS matrix form is consistent with that of the ε rT matrix form as presented in Equation (8): The relationship between ε rT and ε rS is shown in Equation (9): When the piezoelectric stack consists of n basic PMs, and the electric field direction is the same as the polarization direction which is shown in Figure 5 [21], the d 33 mode longitudinal deformation in domain 2.1 is expressed as: The relationship between εrT and εrS is shown in Equation (9): When the piezoelectric stack consists of n basic PMs, and the electric field direction is the same as the polarization direction which is shown in Figure 5 [21], the d33 mode longitudinal deformation in domain 2.1 is expressed as: If the electric field direction is perpendicular to the polarization direction which is seen in Figure 6 [24], the d15 mode shear deformation in domain 2.2 is expressed as: If the electric field direction is perpendicular to the polarization direction which is seen in Figure 6 [21], the d 15 mode shear deformation in domain 2.2 is expressed as: Materials 2020, 13, x FOR PEER REVIEW 8 of 19 Figure 6. The principle of shear deformation caused by the d15 mode.
Assuming that the deformation degree of a basic PM unit is large enough, according to Equations (10) and (11), the deformation effect of n pieces of PM that are loaded with V volts is equivalent to that of a block of PM that is loaded with nV volts.
Supposing that the PIM simulation model is established in a global coordinate system (X, Y, Z), in general, the polarization direction of the PM defaults to the z-axis, so the d33 mode deformation in domain 2.1 is realized by applying an electric field in the z-axis, as presented in Figure 7a. But for the d15 mode deformation in domain 2.2, a rotation coordinate system (RX, RY, RZ) is configured to domain 2.2, making the polarization direction perpendicular to the electric field direction. The definition of a Euler angle α, β, γ is shown in Figure 7b. Equation (11) shows the relationship between the two coordinate systems: Assuming that the deformation degree of a basic PM unit is large enough, according to Equations (10) and (11), the deformation effect of n pieces of PM that are loaded with V volts is equivalent to that of a block of PM that is loaded with nV volts.
Supposing that the PIM simulation model is established in a global coordinate system (X, Y, Z), in general, the polarization direction of the PM defaults to the z-axis, so the d 33 mode deformation in domain 2.1 is realized by applying an electric field in the z-axis, as presented in Figure 7a. But for the d 15 mode deformation in domain 2.2, a rotation coordinate system (RX, RY, RZ) is configured to domain 2.2, making the polarization direction perpendicular to the electric field direction. The definition of the Euler angle α, β, γ is shown in Figure 7b. Equation (12) shows the relationship between the two coordinate systems: According to the polarization direction of domain 2.2 in Figure 7a, the three Euler angles are: In summary, the form of all coefficient matrices in two constitutive equations must be in strict conformity with the above requirements to ensure the model convergence and correct deformation mode. Besides, a rotation coordinate system must be configured in the domain 2.2 to change the polarization direction.

Nonlinear Contact Friction in Solid Mechanics
The contact friction among the four actuation feet and the driving shaft is analyzed by a Coulomb model. Consequently, the tangential traction is obtained from Equation (14), namely, PIM driving force.
where Tt is tangential traction; f is sliding friction; μ is the coefficient of friction; Tn is normal force; Tcohe is cohesion sliding resistance; Tt.max is maximum tangential traction; Tt.crit is critical friction force, Tt.trial is the trial tangential traction; ρt is penalty terms used to express the relationship between the tangential contact force and the penetrations along the tangential direction, it is equivalent to the tangential stiffness; gt is penetration along the tangential direction, which is equivalent to the tangential volume deformation. ρn are the penalty terms used to express the relationship between the normal contact force and the penetrations along the normal direction; it is equivalent to the normal stiffness; gn is the penetration along the normal direction which is equivalent to the normal volume deformation. Moreover, the model nonlinearity is reflected in: the contact interface is According to the polarization direction of domain 2.2 in Figure 7a, the three Euler angles are: In summary, the form of all coefficient matrices in two constitutive equations must be in strict conformity with the above requirements to ensure the model convergence and correct deformation mode. Besides, a rotation coordinate system must be configured in the domain 2.2 to change the polarization direction.

Nonlinear Contact Friction in Solid Mechanics
The contact friction among the four actuation feet and the driving shaft is analyzed by a Coulomb model. Consequently, the tangential traction is obtained from Equation (14), namely, PIM driving force.
where T t is the tangential traction; f is the sliding friction; µ is the coefficient of friction; T n is the normal force; T cohe is the cohesion sliding resistance; T t.max is the maximum tangential traction; T t.crit is the critical friction force, T t.trial is the trial tangential traction; ρ t is the penalty term used to express the relationship between the tangential contact force and the penetrations along the tangential direction, it is equivalent to the tangential stiffness; g t is the penetration along the tangential direction, which is equivalent to the tangential volume deformation. ρ n is the penalty term used to express the relationship between the normal contact force and the penetrations along the normal direction; it is equivalent to the normal stiffness; g n is the penetration along the normal direction which is equivalent to the normal volume deformation. Moreover, the model nonlinearity is reflected in: the contact interface is unknown in advance, including the area size, mutual position and contact state of the contact surface. They change with time (ρ t , g t , ρ n , g n ), which needs to be determined in the solution process. In addition, the contact condition is a unilateral inequality constraint.
The increase of nonlinear components makes the convergence of the model very difficult. It is possible to solve the convergence problem and determine contact traction and friction force by setting the maximum number of iterations in the model solver.

Electrostatics in Piezoelectric Materials
The electrostatic physical field is applied to domains 2.1 and 2.2. The entire physical process can be described by Equation (15): where V is the input voltage signal, P is the polarization vector field, ρ is the space charge density, and ε r is the relative permittivity matrix; it could be ε rT or ε rS , depending on the form of the piezoelectric constitutive equation used. Both E and D can be solved by Equation (15a) and (15b): they are also the input information needed by the piezoelectric constitutive Equations (1) and (5) to solve the strain and stress. Therefore, the main function of electrostatics is to provide the excitation input in the system model.

Setting of Material, Excitation and Mesh
In this section, the material properties of domain 1, domains 2.1 and 2.2, and domain 3 are introduced, respectively. Next, the input excitation is set to achieve the PIM movement. Finally, the model mesh is described.

Setting of Key Piezoelectric Materials and Others
Based on the above theory analysis, the coefficient characteristics matrices of PIC255 in domain 2.1 and domain 2.2 are the key information that must be determined, for example, S E , d, ε rT in the strain-charge equation and c E , e, ε rS in the stress-charge equation.
The coefficient matrices, S E , d, ε rT in the strain-charge are shown in Table 1. The coefficient matrices, c E , e, ε rS in the stress-charge equation are shown in Table 2. The coefficient matrices form are determined by the crystal structure of the PIC255 and polarization. The incorrect coefficient matrix form will cause the model to not be solved, even if the model is successfully solved under certain circumstances, the deformation characteristics of the PM may be incorrect. As for the driving shaft in the domain 1 and fixed axle in the domain 3, C1045 is the driving shaft material according to Reference [24]. Steel AISI 4340 is applied to the fixed axle, and the high rigidity can meet the requirements. Table 3 shows the basic characteristics of the three materials, including PIC255.

Setting of Input Excitation
The laws of motion of the PIM can be obtained from Reference [12] and summarized below: 1.
The movement state of the non-adjacent actuation legs is the same; 2.
For the adjacent actuation legs, when one leg has completed all movements in one cycle, the other leg starts to follow same movements.  As for the driving shaft in the domain 1 and fixed axle in the domain 3, C1045 is the driving shaft material according to Reference [24]. Steel AISI 4340 is applied to the fixed axle, and the high rigidity can meet the requirements. Table 3 shows the basic characteristics of the three materials, including PIC255.

Setting of Input Excitation
The laws of motion of the PIM can be obtained from Reference [12] and summarized below: 1. The movement state of the non-adjacent actuation legs is the same; 2. For the adjacent actuation legs, when one leg has completed all movements in one cycle, the one other leg starts to follow same movements. Figure 8 shows   The four types of excitation functions are listed below: where m 0 and m 1 represent the magnification factor of the longitudinal and shear deformation respectively, depending on number of layers of piezoelectric stack. There, these two magnification factors correspond to n in Equations (10) and (11). The m 0 and m 1 of P-123.01, P-143.01, P-153.01 are presented in Table 4. U is the initial voltage, corresponding to V in the Equation (15). T is the time that requires to complete all movement gesture of one actuation leg. The excitation function needs to be smooth to prevent the model from non-convergence when the magnification factors or initial voltage is set large. If the smooth effect is not good, the smooth transition area needs to be enlarged. According to the official recommended voltage of PIC255, U is set to 250 volts. Moreover, T is one second, the solver time step should be set depending on the model convergence, and the recommended step size is 1 × 10 −2 to 1 × 10 −5 .

Setting of Model Mesh
In response to coupling and nonlinear problems, a fine mesh was utilized to divide the model, but it can bring about a sharp increase in the degree of solution freedom. The computer memory and solution time were greatly increased which seriously reduced the efficiency of model debugging. Therefore, the principle of mesh generation is to minimize the number of meshes without affecting model solution. Figure 9 shows the mesh of the model without driving shaft and the full model. The free tetrahedral mesh was used to partition the two model with super refinement level mesh. In Figure 9a, A is the length of the actuation leg, B and L are the width and height of leg separately. The model with driving shaft is shown in Figure 9b

Results and Discussion
The realization of the PIM motion gesture determines the feasibility of the multi-physics field simulation model of PIM. The motion gesture (d33 mode and d15 mode), driving force, and travel of the PIM were the key factors in the PIM design. All of the above were closely related to the material properties and size (A, B, L) of the PIC255.

Results and Discussion
The realization of the PIM motion gesture determined the feasibility of the multi-physics field simulation model of PIM. The motion gesture (d 33 mode and d 15 mode), driving force, and travel of the PIM were the key factors in the PIM design. All of the above were closely related to the material properties and size (A, B, L) of the PIC255.

Gesture Analysis of the Four Actuation Legs
The longitudinal and shear deformation that was caused by the d 33 and d 15 mode were crucial to the realization of the inchworm gesture. Figure 10 shows the movement gesture change of the actuation legs in the inchworm model without a driving shaft. The gesture changes of four actuation legs are presented in Figure 10a-d; in addition, the leg gesture at 0 s was exactly the same as the gesture at 1 s. It can be concluded that the movements of non-adjacent actuation legs are the same at each time point. For a more detailed analysis, one-leg gesture change is highlighted in Figure 10e, as for adjacent actuation legs, when the "walk" action of one actuation ends, the other one starts to "walk".

Driving Force Analysis
The driving force was equal to the contact tangential traction. If the driving force is greater than the maximum shear load (MSL), the piezoelectric stack will be broken in the shear direction. The maximum driving force was equal to the maximum shear load.

Driving Force Analysis
The driving force was equal to the contact tangential traction. If the driving force is greater than the maximum shear load (MSL), the piezoelectric stack will be broken in the shear direction. The maximum driving force was equal to the maximum shear load. Table 5 presents the standard size and MSL of P-123.01, P-143.01, and P-153.01 which were obtained from PI Motion.

Driving Force Analysis
The driving force was equal to the contact tangential traction. If the driving force is greater than the maximum shear load (MSL), the piezoelectric stack will be broken in the shear direction. The maximum driving force was equal to the maximum shear load. Table 5 presents the standard size and MSL of P-123.01, P-143.01, and P-153.01 which were obtained from PI Motion. Equation (17) is the absolute error formula. Figure 11 depicts the contact tangential traction of the above three products over time in the simulation, the MSL of P-123.01, P-143.01, and P-153.01 were 40.611N, 203.286N, and 291.657N respectively, the corresponding absolute error were 1.50%, 1.65%, and 2.78% and which are shown in Table 5. Therefore, the maximum shear load increases with the increase of PM size; this conclusion is also consistent with the MSL, as shown in Table 5

Travel Analysis
The travel was the displacement of the driving shaft in the shear direction during the action completed period (T) as illustrated in Figure 12. The travel data from manufacturer and simulation are compared in Figure 13, P-123.01′s travel in simulation was 1.6μm, and there was an error of 0.6 μm compared with the manufacturer's standard data. The travel of P-143.01 and P-153.01 were 0.9 μm and 3.02 μm, and the errors were 0.1 μm and 0.02 μm, respectively. In addition, a large number of meshes were-needed for high-solution accuracy. Model calculations with many meshes and coupling of multiple physical fields has high requirements for computer performance and needs a long solution time.

Travel Analysis
The travel is the displacement of the driving shaft in the shear direction during the action completed period (T) as illustrated in Figure 12. The travel data from manufacturer and simulation are compared in Figure 13, P-123.01 s travel in simulation was 1.6µm, and there was an error of 0.6 µm compared with the manufacturer's standard data. The travel of P-143.01 and P-153.01 were 0.9 µm and 3.02 µm, and the errors were 0.1 µm and 0.02 µm, respectively. In addition, a large number of meshes were needed for high-solution accuracy. Model calculations with many meshes and the coupling of multiple physical fields have high requirements for computer performance and need a long solution time.

Conclusions
In this paper, a multi-physics coupling model of an inchworm motor was established based on PI Motion's advanced inchworm motor products. The motion gesture and performance parameters of the inchworm motor were the research priority for the motor design and performance estimation. Three different sizes of products, including P123.01, P143.01, P153.01, were simulated with 250 V initial voltage in one second. The key piezoelectric material was PIC255. Two conclusions are as follows: The d33 and d15 mode deformations of PIC255 were perfectly realized by simulation. That is, the inchworm's gesture characteristic was perfectly realized, and it is highly consistent with PI Motion's products. The realization of motion gesture proves that the multi physical field coupling simulation is feasible for assisting the PIM design.
The maximum driving force error of P123.01 between the simulation data and the official data was 2.8%, the errors for P143.01 and P153.01 were 1.50% and 1.65%, respectively. The driving forces of the three types of products obtained from the simulation were close to those provided by the actual official data. In addition, the driving force increases with the increasing piezoelectric material size. As for the PIM travel, there was a maximum travel error of 0.6 μm between the simulation and the official standard for the P123.01, the errors for P143.01 and P153.01 were 0.1 μm and 0.02 μm. The travel data from the simulation and the official standard was well matched. So, the performance of the PIM can be directly obtained by multi-physics field coupling simulation.
In conclusion, the multi-physics field coupling simulation model involving solid mechanics and electrostatics is suitable for PIM performance evaluation and assisting motor design.
Author Contributions: B.Z. contributed to the conception of the study and performed the calculations; R.F. performed the simulation and drafted and revised the manuscript; W.S. analyzed, reviewed, and edited manuscript. All authors have read and agreed to the published version of the manuscript.

Conclusions
In this paper, a multi-physics coupling model of an inchworm motor was established based on PI Motion's advanced inchworm motor products. The motion gesture and performance parameters of the inchworm motor were the research priority for the motor design and performance estimation. Three different sizes of products, including P123.01, P143.01, P153.01, were simulated with 250 V initial voltage in one second. The key piezoelectric material was PIC255. Two conclusions are as follows: The d 33 and d 15 mode deformations of PIC255 were perfectly realized by simulation. That is, the inchworm's gesture characteristic was perfectly realized, and it is highly consistent with PI Motion's products. The realization of motion gesture proves that the multi physical field coupling simulation is feasible for assisting the PIM design.
The maximum driving force error of P123.01 between the simulation data and the official data was 2.8%, the errors for P143.01 and P153.01 were 1.50% and 1.65%, respectively. The driving forces of the three types of products obtained from the simulation were close to those provided by the actual official data. In addition, the driving force increases with the increasing piezoelectric material size. As for the PIM travel, there was a maximum travel error of 0.6 µm between the simulation and the official standard for the P123.01, the errors for P143.01 and P153.01 were 0.1 µm and 0.02 µm. The travel data from the simulation and the official standard was well matched. So, the performance of the PIM can be directly obtained by multi-physics field coupling simulation.
In conclusion, the multi-physics field coupling simulation model involving solid mechanics and electrostatics is suitable for PIM performance evaluation and assisting motor design.