Mechanical Response of MEMS Suspended Inductors under Shock Using the Transfer Matrix Method

MEMS suspended inductors are susceptible to deformation under external forces, which can lead to the degradation of their electrical properties. The mechanical response of the inductor to a shock load is usually solved by a numerical method, such as the finite element method (FEM). In this paper, the transfer matrix method of linear multibody system (MSTMM) is used to solve the problem. The natural frequencies and mode shapes of the system are obtained first, then the dynamic response by modal superposition. The time and position of the maximum displacement response and the maximum Von Mises stress are determined theoretically and independently of the shock. Furthermore, the effects of shock amplitude and frequency on the response are discussed. These MSTMM results agree well with those determined using the FEM. We achieved an accurate analysis of the mechanical behaviors of the MEMS inductor under shock load.


Introduction
The development of MEMS technologies introduces new approaches and designs to the fabrication of inductors in order to attain higher quality factors that are required for typical microwave circuits [1]. A MEMS inductor can be separated from the substrate by etching the sacrificial layer or the substrate, which is called the MEMS suspended inductor [2][3][4][5][6]. Although this structure decreases substrate loss in terms of electrical performance [1], it is sensitive to external forces, as the spiral is movable in terms of structural dynamics, just as MEMS capacitive accelerometers are [7]. The most severe external force that MEMS devices are subjected to during operations is mechanical shock. In particular, the normal shock acceleration during the drop process is around 500 g (g is the acceleration of Earth's gravity); during the launch step of a high-speed vehicle, the devices are subjected to an acceleration that is much greater than 10,000 g [7,8]. This means that in some probable working situations, under such high-g shock loads, excessive deformation of the inductor might result in the degradation of the electrical properties, with some parts yielding or even fracturing, which is the most serious failure. There have been some research studies on MEMS suspended inductors in shock environments. Ribas et al. [2], Lin et al. [5] Hsieh et al. [6], and Dahlmann et al. [9], respectively, used the finite element program called ANSYS to examine the deformation of MEMS suspended inductors under shock acceleration, in addition to analyzing the electrical properties. In these studies, as the shock loads were below the drop level, the deformations of the inductors were almost negligible, and the researchers did not investigate further. The finite element method (FEM) mentioned above is generally used when the problems are too complicated to be solved by an analytical method, such as the deformation of three-axis MEMS gyroscopes [10].
On the other hand, researchers have proposed some analytical models of MEMS devices under shock [11]. Jiang et al. [12], Heish et al. [6], and Dahlmann et al. [9] compared the inductor to a cantilever beam and considered that the maximum displacement of the free end of the beam was equivalent to the maximum deformation of the inductor spiral. These treatments were rough and their results were incorrect. Likewise, Srikar et al. [13] proposed another equivalent model for treating a MEMS device as a single-degree-of-freedom (SDOF) system and obtained damage criteria formulations based on displacement and stress. Xu et al. [14] and Peng et al. [15] have each proposed more accurate formulations for describing shock acceleration based on the SDOF model. Wang et al. [16], Lian et al. [17], and Fathalilou et al. [18] have analyzed the shock resistance performance of MEMS gyroscopes and MEMS shock switches using a dual-vibrator model, an extension of the SDOF model, respectively. Fang et al. [19] and Singh et al. [20] combined both the SDOF model and the microbeam model to analyze the shock resistance of MEMS micromirrors and MEMS switches. Sundaram et al. [21] proposed a combined experimental-analytical approach by focusing on equating key parts of the device to the SDOF systems and they investigated the failure of a tunable diffraction grating under shock and vibration. Xu and Li [22][23][24] have carried out many meaningful studies. They equated the inductor to an SDOF system to obtain an equivalent acceleration, considered the inductor as a bar system, determined its deformation and stress by solving the super-stationary equations, and analyzed the mechanical reliability of the structure. The natural frequency of the structure, however, was established by the FEM; therefore, it was not a fully analytical solution. Moreover, they did not discuss the relationship between the shock and the response in detail, but only qualitatively. In fact, most analytical studies on microstructure mechanics have focused on a class of devices represented by microbeams [13,[25][26][27][28], with little study on MEMS suspended inductors. Therefore, the situation calls for an accurate analytical method to analyze the response of the MEMS inductor under shock and to determine whether the inductor has failed since understanding the dynamics of a device is critical to determine its failure.
A MEMS suspended inductor could be regarded as a multibody system composed of rigid and flexible bodies in structural mechanics. In the research on multibody system dynamics, Rui et al. [29][30][31] systematically proposed the transfer matrix method for linear multibody systems (MSTMM) after more than 20 years of continuous refinement. This method introduces a new approach to analyzing the dynamics of a multibody system by establishing transfer relationships between the state vectors of each component. In addition, the authors applied the method to the dynamic analysis of laser gyros, piezoelectric actuators, rocket boosters, etc [32][33][34][35][36][37].
This paper describes in detail the use of the MSTMM to investigate the dynamic response of the MEMS suspended inductor under high-g shock. The natural frequencies and mode shapes were obtained according to the geometry and material parameters, and the dynamic response of the structure was solved by modal superposition to obtain the deformation and stress, as well as the maximum shock acceleration when the yield strength was reached.

Simplification of the Structure of the Inductor
Before the mechanical analysis, the real structure should be simplified, based on the characteristics of the individual problem, in order to skip over the unnecessary details and retain the basic features. The following describes the assumptions and simplifications used in building the mechanical model of the MEMS suspended planar inductor. A typical MEMS suspended planar inductor is made up of three parts: a spiral, a substrate, and two pillars, as shown in Figure 1. The spiral could be square-polygonal or circular in shape. In this paper, we investigate a 1.5-turn square spiral made of copper with a silicon substrate. The spiral was supported by two pillars on the substrate, which also connected the signal transmission lines. MEMS inductors, which are designed in the same layout but manufactured using different processes, might not be identical in several aspects, including metal thickness, surface flatness, silicon resistivity, silicon dioxide thickness, residual stress [38], and so on. Here, it is assumed that the three-dimensional dimensions of the metal spiral are the correct sizes in the design layout, in the sense that the spiral strips are cuboids in geometry. Considering the shock, the load causes a sudden acceleration that acts on the object, which can be described by an acceleration-time curve and has three basic characteristics, e.g., strength, shape, and frequency. Shock loads in different environments are generally highly varied, even when generated by impact equipment used in the laboratory. The shock acceleration can usually be approximated by a half-sine pulse of the shock peak, which allows for the uniform analysis of a large class of shock problems, as shown in Figure 2. For MEMS devices, regardless of how the shock load is generated, it can be thought of as a distributed acceleration pulse acting on the package and substrate, characterized by a short duration (the package is assumed to have minimal influence on the strength and shape of the shock load [13]). Furthermore, because the damping of the system does not absorb enough energy in the duration, the damping effect could be ignored. In the following, the MSTMM is used to analyze the problem of the dynamic response of the MEMS inductor under the shock load. The analysis could be organized into several steps. First, the system is decomposed into six components; next, the transfer equation of each component can be established. Each strip satisfies the definition of a bar so that it can be reduced to its axis. The connected parts of the adjoining strips are represented by the intersection of the axes, and the action points of the loads are transferred to the axes, as illustrated in Figure 3. The support pillars that connect the spiral to the substrate cannot be equated to its axis. In terms of this problem, it is considered that the pillars, together with the substrate, impose fixed constraints on the spiral due to their minimal deformation under shock loads, thus allowing them to be treated as rigid bodies.
It is clear that the spiral is the movable part of the structure and the direction perpendicular to the plane (in which the spiral is located) is the most prominently affected by the shock load. When subjected to the load in the normal direction, each bar undergoes deformation, which is mainly a combination of bending and torsion. Even though stress and deformation are more complicated, as long as the deformation is small and the bar remains linearly elastic, it can be assumed that the various basic deformation forms are independent of each other, which is the principle of the superposition of the combined deformation of the bar.

Free Vibration Characteristics of the Inductor
The next step is to determine the vibration characteristics of the spiral, including free and forced vibration characteristics. The free vibration characteristic is analyzed first; this includes the natural frequencies and mode shapes. Consequently, the transfer matrix of each component needs to be obtained by the MSTMM. Moreover, considering that the six bars have the same physical and material parameters, except the length, only the transfer matrix of the first bar needs to be derived, and the others can be acquired in the same manner.
The transfer equation of the bar is derived from the motion-governing equations of bending and torsion. Equation (1) is the partial differential equation, where the bar is bent, which describes the lateral free vibration of the bar, and is called the Euler-Bernoulli beam equation. In addition, according to the MSTMM, each element has its own stationary inertial Cartesian coordinate system. The origin is located at the endpoint that is first encountered along the transfer direction [31].
wherem is the linear density of the bar and EI is the bending stiffness of the beam; it is more convenient to set z(x, t) = Z(x)e iωt Substituting the above into Equation (1), in the modal coordinates, there is The general solution of Equation (2) is where β = 4 mω 2 (EI) . Considering the knowledge of geometry and the mechanics of materials, where Θ y denotes the angle displacement and Q z and M y denote the interior force and torque in the modal coordinates. (4), and arranging the outcome into a matrix,

Substituting Equation (3) into Equation
the above can be abbreviated as where Z(x) denotes the state vector at x in the modal coordinates, C is the matrix of the integral constants. The state vector Z is a column vector representing the mechanical state of any point in a multibody system; the elements are the (angular) displacement, the internal force, and the internal moment at that point.
In order to formally eliminate the integral constant vector C, the operations can be done as follows. when In particular, when x = l (l is the length of the continuum axis of a particular component), Z(l) = U(l)Z(0). More generally, it can be written as where U is the transfer matrix of the bending bar, which does lateral free vibration. The subscripts O and I indicate the output and input points of the component, respectively.
For the bar bending problem, there is Equation (9) represents the partial differential equation that arises when a round bar is twisted.
where ρ is the volume density, G is the shear modulus, and J p is the polar moment of inertia. In this problem, some changes must be made to the parameters in Equation (9) because the cross-section of the bar is not round but square. According to the related research of elasticity, J p in the torsional stiffness GJ p is replaced by the equivalent polar moment of inertia J s of a rectangular cross-section bar. Moreover, J p in ρJ p is replaced by the polar moment of inertia J o of the rectangular cross-section. Similarly, in the modal coordinates, The general solution to the above is where γ = ω ρJ o GJ s .
Considering the mechanical relationship of materials, Substituting Equation (11) into Equation (12) and arranging the outcome into a matrix, Eliminating the integral constants matrix and writing it in the form of Equation (7), By combining Equations (8) and (14), we can derive the transfer matrix equation for a bending and twisting bar as follows: where In fact, the transfer equations of all types of elements can be written in this form. For convenience, the endpoints of all components are marked counterclockwise from 0 to 6, and the corresponding state vectors are marked Z 0 to Z 6 , so that where U j denotes the transfer matrix of the j-th component and Z j denotes the state vector of the j-th point. Note that the output point of one component is exactly the input point of the next component, which means that the state vectors of these two points should be equal, taking into account the rotation of the coordinate system. For a chain system, such as the inductor spiral, the total system transfer equation can be obtained by assembling all transfer equations in a counter-clockwise order. Thus, the total transfer equation of the spiral is where U all is the total transfer matrix of the system and is a 6 by 6 matrix, H represents the direction cosine matrix that is rotated 90 degrees counterclockwise. It can be seen that the total transfer equation of the system only involves the state vectors of the boundary points and not the state vectors of the connection points. Moreover, the elements in the total transfer matrix U all are simply functions of the structure parameters and the natural frequency. By substituting the boundary conditions into Equation (17), the homogeneous linear equations, which are called characteristic equations, can be extracted. In general, half of the elements in the boundary point state vector are unknown and half are known. In the study, when both ends are fixed, Substituting the above into Equation (17), the corresponding elements in Z 0 , Z 6 , U all form the characteristic equation.
The force Q z and the torque M x , M y must exist, So the determinant of the coefficient matrix U should be zero, detU = 0. The natural frequency ω k (k = 1, 2, · · · ) can be obtained by solving the determinant in which the only unknown quantity is the natural frequency of the system. In the case that Q z,0 = 1, the values of M x,0 , M y,0 can be derived by solving non-homogeneous linear equations, which are part of the elements of Equation (19). In this way, the state vector of the start point, Z 0 , is obtained, and the values of Z 1 to Z 6 can also be obtained by substituting Z 0 with Equation (16). The state vector at any position of each element is given by Z j (x) = U j (x)Z j−1 . Moreover, the displacements and angular displacements in all of the Z j (x) constitute the modal shapes of the system.

Dynamics Equation of the Inductor
The following describes how to solve the response of the spiral under shock acceleration. The dynamic equations of each component can be written in the following form [31]: where M j and K j are the parameter matrices, which can be called the mass matrix and stiffness matrix, following the habits of vibration mechanics. v j is the column matrix of displacement (angular displacement). v j,tt and v j,t denote the two-order derivative and one-order derivative of v j , with respect to time t. f j is the column matrix of the force (torque) acting on the j-th body (j is the body index in the multibody system). Likewise, air damping is not taken into account. The body dynamics equation of the bending and twisting bar can be obtained by arranging the partial differential Equations (1) and (9), taking into account the force to the matrix.
where a is the shock acceleration. The total body dynamic equation of the system can be easily obtained by arranging the body equations of all components in turn, where

Orthogonality of Eigenvectors
In the case of free vibration, the elements in the column matrix of displacement v can be represented in modal coordinates, v = Ve iωt . Then V , v in the modal coordinates can be obtained, as shown in the following formula, which is denoted as the eigenvector of the system.
where k is the sequence number of the natural frequency of the structure. According to the characteristics of the parameter matrices M and K, it can be proved that the two pairs of inner products are equal, respectively, which is called the symmetry of eigenvectors, with respect to M and K [31].
The sufficient mathematical proof of the symmetry of the eigenvectors is given in Appendix A.
For the free vibration of the k-th mode, v = V k e iω k t . Substituting v into Equation (22), in the case of free vibration, so Eventually, there is when ω k 1 = ω k 2 , MV k 1 , V k 2 = 0. The same procedure can be easily adapted to prove that KV k 1 , V k 2 = 0(ω k 1 = ω k 2 ). These show that the mode shapes with respect to different frequencies are orthogonal to each other for the mass and stiffness matrices. For convenience, they can be written as where

The Dynamic Response of the Inductor
Having established the orthogonality of the eigenvectors V , the displacement response v can be expressed by the product of the generalized coordinate q(t) and the eigenvectors V , as seen below, according to the modal superposition method [31].
Substituting Equation (29) into the total body dynamics Equation (22), then Taking the inner product on both sides with the eigenvector V k and using its orthogonality, which is a single-valued function that only depends on time t. Using Duhamel's integral to obtain the solution of Equation (31), on the condition that the initial position and initial velocity of the system in this problem are both zero, it holds that Finally, the dynamic response of the system is

The Maximum Z-Axis and Angular Displacements
After obtaining the dynamic response of the inductor, the maximum displacement problem and von Mises stress are discussed. Taking out the Z-axis displacement function z j (x) of each element from the response matrix v, there is To calculate f k (t) in the above equation, we substitute the function expression of the half-sine wave shock pulse, a(t) = a 0 sin ω f t(t ≤ τ), into Equation (32), where N k depends on ω k . Substituting the above into Equation (35), there is It can be seen that, for any point on the spiral, its response is similar to the forced vibration of an undamped single-degree-of-freedom system. Calculating the integral, each component in Equation (37) can be expressed as Likewise, the angular displacement θ x,j of each element has a similar mathematical form.
It is generally accepted that the maximum response caused by the shock load is more meaningful than other responses throughout the entire reaction process. It is the right time to discuss where and when the maximum displacement and stress occur. The displacement response is the sum of the contribution of each mode shape. In most cases, the lowest frequency portion makes the largest contribution, and the higher frequency portions tend to decrease significantly [39]. Moreover, as can be seen from Equations (38) and (41), the roles, in terms of the functions of position and time in the contribution of the modal shapes, are independent. Therefore, it is only necessary to find the position and time corresponding to the maximum of the first-order portion. This specific position can be obtained by plotting the position function, and the specific moment needs to be discussed on a case-by-case basis. First, we discuss the case of t ≤ τ in Equation (39), take the derivative of the function, and set it equal to zero; we can obtain cos ω f t = cos ω k t. Considering that n is a positive integer); in the second case, the problem is simpler, where t should meet t = π 2 ( 1 ω f + 1 ω k ).

The Maximum von Mises Stress
With the column matrix displacement v obtained, the moment and stress can be calculated according to the relationship between displacement, moment, and stress. For each bar component, the bending moment on the cross-section at position x is The maximum tensile stress appears at the long side of the cross-section.
Moreover, the torque on the cross-section at position x is The maximum shear stress appears at the midpoint of the long side of the cross-section.
Furthermore, the maximum von Mises stress appears at the midpoint of the long side of the cross-section.
Substituting Equation (37) with (41) into the above equation, and making the necessary arrangements to obtain where b and h are the lengths of the long and short sides of the beam cross-section, respectively.
The position and time of the maximum von Mises stress are determined by the position and time of the maximum contribution of the first-order modal shape. So, only considering the case of n = 1, Equation (47) becomes This formula demonstrates that the maximum von Mises stress occurs at the same time as the displacement, and the position is determined by graphing the function.

Results and Discussion
This section mainly includes the results obtained by the MSTMM and the FEM to deal with the response of the MEMS suspended inductor under shock. The natural frequencies of the inductor were calculated, as well as the displacements and von Mises stresses under different loads. The accuracies of these data were verified by FEM. The material parameters used are listed below, including the parameters involved in the previous section. The Young's modulus of copper is 128 GPa, the Poisson's ratio is 0.34, and the density is 8900 kg/m 3 . The cross-section of the strips had a width(b) of 20 µm and a height(h) of 10 µm. The spiral was 1.5 turns with a maximum outer diameter of 200 µm and a gap of 10 µm.

Natural Frequencies of the Inductor
The natural frequencies obtained by FEM pertained to all degrees of freedom in the system. For this problem, only the frequencies corresponding to the z-direction mode shapes were required. Therefore, these frequencies were selected and compared with the results obtained by MSTMM. As can be seen from Table 1 below, the natural frequencies of the MEMS inductor were relatively high, partly because the mass was generally very small. The errors between the data calculated by the two methods were small and should be due to the simplifications of the model.

Responses of the Inductor Under Shock
A half-sine pulse load was selected with an amplitude of 10,000 g, a frequency of 10 kHz, and a direction along the negative z-axis, which meant that the shock load lasted 50 µs. Figure 4a,b shows, respectively, the displacement and the maximum Von Mises stress of each cross-section of the strips at a given moment (actually, these moments were the times at which the maximum values occurred). The responses of the displacements obtained by both methods were in good agreement. It is apparent that the maximum response of vertical displacement appeared at the intersection of components B3 and B4, namely P3. Moreover, the displacement results for B3 were larger than the FEM results (maximum relative error is 7.5%), which should be due to the approximation of the torsion equation.
For the von Mises stress, however, the situation became a little more complicated. Excluding all intersections and the surrounding region, the stresses exhibited relatively good agreement. However, at the intersections, the results obtained by FEM became particularly small. It indicated that stress concentrations appeared at these corners. The variation in the cross-section caused the stresses to be concentrated in the regions near the corners, which meant that the maximum values of the stresses were not at the midpoints of the long sides of the cross-sections. The previous method of calculating the stress was no longer valid. To solve this problem, a method was proposed to estimate the concentrated stress by vector-summing the shear and normal stresses of the two cross-sections of the two adjoining components corresponding to the corner point, and then calculating the von Mises stress; the results are shown in Table 2. It can be seen that the relative errors were small at P1, P4, and P5, where the stress concentration effects were more obvious; the relative errors were very large at P2 and P3, where the effects were not so obvious. Moreover, the results obtained by FEM were related to the mesh density: the denser the mesh, the higher the stress values. Thus, the maximum von Mises stress appeared at the intersection of elements B4 and B5, namely P5. The stresses at the first fixing points were very high as well.   Next, we present the results relating to the time at which the maximum value occurred. Here, an additional shock with a frequency of 70 kHz was chosen to verify the correctness of the MSTMM and to determine the timings of the maximum responses at different frequencies. Figures 5 and 6 present the time histories of the vertical displacements at position P3 and the von Mises stress at corner P5 on the spiral when subjected to shocks with frequencies of 10 kHz and 70 kHz, and they indicate that the shapes of the time histories of the responses were almost identical at all points if the response at each point was taken as the absolute value and normalized, which was exactly the function curve of |φ(t)| given in Equation (39). This also meant that the maximum displacement occurred almost at the same time as the maximum stress. Looking first at Figure 5, the maximum response occurred during the forced vibration phase when the frequency of the shock was lower than the first-order natural frequency. The time of the maximum response obtained by the MSTMM coincided with the first peak of the curve, whereas the time obtained by the FEM corresponded to the second peak. This was due to the slight difference in the natural frequencies obtained by the two methods, which in turn affected the n in the expression of the times of the maximum being 1 and 2, respectively. The relative error of the maximum displacement response was 7.6% and the relative error of the stress was 4.0%. Looking at Figure 6, the maximum response occurred in the free vibration phase due to the ratio of the force and the first-order natural frequency being greater than 1. The relative error of the maximum displacement response was 5.1% and the relative error of the stress was 5.0%. In this subsection, we determine the time and position at which the maximum response of the inductor under shock load occurs, which is very important for the following analysis.    Figure 8 shows the maximum displacement and maximum von Mises stress of the inductor under half-sine shocks with the same amplitude of 10,000 g but different frequencies.

Effect of Shock Frequency on Responses
Similarly, it can be seen that the displacement and the stress curves exhibit nearly identical shapes when normalized. The maximum displacement and stress both occurred at a frequency of 32 kHz. This means that the most dangerous situation would appear when the frequency of the force is about 32 kHz, which corresponds to a duration of 15.6 µs; in this case, an amplitude of 36,000 g would cause the von Mises stress to exceed the yield strength of the body copper by 200 MPa. It is also important to note that the displacement must be less than the suspended height. One more pillar could be added at P3 to prevent the spiral from fracturing. The advantage of the analytical method was shown when dealing with problems of different external forces, where only a few parameters related to the external forces needed to be changed. In contrast, the FEM had to repeat the entire process. For the cases in this subsection, the MSTMM took only 5 min, while the FEM took over 30 min.
From the previous discussion, it can be confirmed that it is possible to calculate the maximum displacement and maximum von Mises stress (within the yield strength) of the MEMS inductor under a shock of any amplitude and frequency. The responses at different frequencies with an amplitude of 10,000 g can be found in Figure 8a,b, respectively, and the values can be obtained by multiplying the corresponding ratios of the amplitudes.

Conclusions
In this paper, the response of the MEMS suspended inductor under high-g shock was investigated using the transfer matrix method of the linear multibody system, MSTMM. The free vibration characteristics of the system and the dynamic response under shock load were obtained. It was theoretically demonstrated that the positions where the maximum displacement and von Mises stress occurred can be determined independently of the shock load. Moreover, the procedure used to find the time of the maximum was explained; the method of identifying the maximum von Mises stress in the case of stress concentration was also discussed. The effect of the shock amplitude on the response was linear, and the response curves of the system were given for different shock frequencies. Therefore, the maximum response value to the shock load at a certain frequency at any amplitude could be quickly determined by looking up values in the frequency-response table and multiplying them by the amplitude ratios. This approach helps identify where the inductor may yield or fracture. Our work could serve as the foundation for designing shock-resistant MEMS suspended inductors and as a mechanical basis for further research on the relationship between changes in inductor electrical properties and shock loads.