# Inverse Eigenvalue Theory-Based Rigid Multibody Modeling Method of Complex Flexible Structures in Large-Scale Mechanical Systems

Authors: Jiulin Wu; Lizhan Zeng; Bin Han; Xin Luo; Xuedong Chen; Wei Jiang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Increasing attention is paid to modeling flexibility of individual components in the multibody simulation of large-scale mechanical systems. Nevertheless, the high model order of common methods such as FEA restricts efficient explorations, especially in dynamic design and iterative optimization. In this paper, a rigid multibody modeling strategy (RMMS) with low DOFs and explicit physical meaning is proposed, which directly discretizes a continuous structure into a number of rigid finite elements (RFEs) connected by spring-damping elements (SDEs). In the RMMS, a new identification method from the perspective of the inverse vibration problem is particularly put forward to resolve the parameters of SDEs, which is crucial to the implementation of RMMS in complex flexible structures. With decoupling and linearization, this nonlinear problem is transformed into solving the incompatible linear equations inRn2 vector space based on vectorization operator and Kronecker product, and optimal parameters are obtained by calculating the Moore–Penrose generalized inverse. Finally, the comparison of the experimental results with the simulated ones by the RMMS strongly validates the feasibility and correctness of the RMMS in predicting the dynamic behaviors while with few DOFs and explicit physical meaning; the application in a lithography system exhibits the applicability of the RMMS for dynamic modeling of large-scale mechanical systems.


---

## Body

## 1. Introduction

Large-scale mechanical systems such as precision equipment, vehicles, aircraft, and robotics, which comprise plenty of components, are commonly modeled as multibody systems (MBSs) [1]. With the increasing demand for simulation accuracy, the requirement to model the flexibility of weak-rigid components in these systems has become much more stringent [2, 3]. Imposing the FEA method in this process serves as the most widely adopted approach, however, yielding a complex model with high DOFs that can easily reach millions or even higher level [4]. The huge computing time and high resource requirement make it inadequate for dynamic analysis and iterative optimization. So, a model with less DOFs which still can catch the nonnegligible dynamic effects caused by structural deformations is necessary.Considerable efforts have been made in model order reduction (MOR) technologies to effectively lower the model order while preserving the behavior of the original system as possible. Representative structural MORs contain (i) condensation methods, including Guyan-Irons static reduction [5], dynamic condensation [6], and improved reduced system [7], which reduce the model by preserving master DOFs while excluding slave ones, (ii) the Krylov-based methods (also known as moment matching) [8], which is essentially a Pade approximation of transfer functions as the Krylov subspace spans derivatives of the transfer functions, (iii) Component Mode Synthesis (CMS) methods, including the fixed-interface method [9], free-interface method [10], and hybrid method [11], in which the reduced subspace is formed by hybrid physical and modal basis vectors, and (iv) other MORs such as Truncated Balanced Realization (TBR) [12] and Proper Orthogonal Decomposition (POD) [13]. These MORs are widely applied but still encounter difficulties since most of them seek to reduce the mathematical models (system matrices) rather than original physical models [14]. Firstly, most of the reduced coordinate spaces have no explicit physical meaning, so that models in reduced space may not be recovered into original physical space, especially during optimization. Secondly, the reduced models have poor reusability and extensibility that a small change in original structure may need a repeated MOR.Rui et al. and Yang et al. [15, 16] proposed the transfer matrix method of the multibody system (MSTMM), which processes the following features: without global system dynamic equations, low model orders, high computational efficiency, and high programming and then developed a virtual design software named MSTMMSim to realize the rapid visual simulation and virtual design for complex multirigid-flexible-body mechanical systems. Kruszewski et al. [17, 18] proposed an alternative way named the rigid finite element method (RFEM) to model the flexibility, which can overcome the abovementioned difficulties of MORs since equivalent rigid MBSs with low DOFs and explicit physical meaning were constructed. The RFEM discretizes a continuous structure into a finite number of RFEs, which reflect the mass (inertial) features of the body, by two division steps. These RFEs are connected by means of massless and nondimensional SDEs, which reflect the elastic and damping features. Based on this concept, Wittbrodt et al. [19, 20] systematically summarized and gave a detailed derivation of the RFEM, and then proposed a modified rigid finite element method (MRFEM), in which the translation DOFs of SDEs are constrained to gain models with more lower order.The RFEMs (include RFEM and MRFEM) have been widely employed in modeling beam-like flexible structures, whose longitudinal dimension is evidently larger than the transverse ones. Wittbrodt et al. [19–21] utilized the RFEMs in the dynamic analysis of A-Frame of ships, telescopic, and offshore mechanical structures such as BOP transportation gantry, offshore crane, and pipeline laying machine. Adamiec-Wójcik et al. [22–25] applied the RFEM in dynamic analysis of plane manipulators, marine ship risers, and flexible ropes to solve the vibration and drive control issues. Krukowski and Maczynski [26] constructed a rigid finite element model of a pedestal offshore crane and analyzed its dynamic response considering flexibility in selected directions. Cao et al. [27] applied the RFEM in the modeling of a ball-bearing rotor system, and comparison between measurement and simulation shows the ability of the model to predict the dynamic behavior of the systems. Li and Jiang [28] built the rigid finite element model of a ball screw system, and experimental results show that the model can be used for vibration analysis and control of the screw system. In addition, the RFEMs were also extended to model 2D flexible structures, such as plates and shells [29], reinforced concrete slabs [30], and collecting electrodes [31]. Additional applications based on the RFEMs can be found in the literature [32–36]. From these studies, it can be seen that although the RFEMs were extensively adopted, they were mainly applied to relatively simple flexible structures. The reason lies on the fact that the application of RFEMs highly depends on the parameter identification of SDEs, which is originally proposed for beam-like structures based on generalized Hooke’s law [19], and then extended for structures like plates and shells. Thus, the RFEMs are inapplicable to model complex flexible structures which possess characteristics such as intricate topology, irregular cross sections, and material inhomogeneity, due to the difficulty in parameter identification of SDEs.In this paper, a novel RMMS is proposed to model complex flexible structures. The RMMS preserves the advantages of RFEMs that have been proven, including (i) simplicity and high numerical efficiency, which are reflected in the description of a system with few DOFs, (ii) modeling of rigid and flexible bodies in a uniform approach, and (iii) applicability for the analysis of both small and large deformations, as well as for static and particularly dynamic analysis [19, 34, 37]. The generated rigid MBSs of complex structures have low model order and explicit physical meaning to effectively predict and optimize the system dynamics. The main challenge in implementation of the RMMS lies on the parameters of SDEs; hence, this paper lays emphasis on this parameter identification problem and proposes a new method to resolve it from the perspective of the inverse eigenvalue problem. This nonlinear problem is eventually decoupled and linearized into solving the incompatible linear equations in Rn2 vector space, and optimal parameters are obtained by calculating the Moore–Penrose generalized inverse. The feasibility and validity of the parameter identification method and the proposed RMMS are finally demonstrated by experiments; and the RMMS is applied in dynamic modeling of a large-scale mechanical system.This paper is organized as follows. Section2 outlines the entire modeling strategy. Division and pre-evaluation, which prepare for the parameter identification, are presented in Section 3. Section 4 elaborates the parameter identification of SDEs, which should be emphasized. In Section 5, the experiments are conducted and the results are compared with the simulated ones. Conclusions are drawn in Section 6.


## 2. Description of the RMMS

For modeling of complex flexible structures, a forward train of thoughts as in RFEMs is hard to be implemented since it is quite difficult to directly identify the parameters of SDEs by common means of theoretical formulas based on material mechanics in ref. [19]. Thus, we propose the RMMS, as described in Figure 1, to model the structural flexibility, which yields an equivalent rigid MBS with few DOFs and explicit physical meaning.Figure 1
A schematic diagram of RMMS: (i) input flexible body, (ii) parametric MBS, (iii) output equivalent rigid MBS, and (iv) FE modal results.In the RMMS, firstly, a flexible structure is directly divided into RFEs instead of two division steps in RFEMs. According to the topological correlations of RFEs in the original structure, SEDs are employed to connect the RFEs. Considering that the parameters of SDEs characterizing spring and damping features are unknown, an equivalent parametric MBS is thereby established. Secondly, FE modal results of the flexible body are calculated, pre-evaluation can be performed, which superiors to RFEMs, and the target eigenpairs are generated in combination with the division scheme. Thirdly, dynamic parameters of SDEs are identified by means of the inverse eigenvalue theory (IET), in which the vectorization operator and Kronecker product are introduced to transform the nonlinear inverse eigenvalue problem into solving incompatible linear equations. Then, the final equivalent rigid MBS is obtained by substituting the identified parameters.It should be emphasized that this parameter identification problem can finally be concluded as an inverse eigenvalue problem on account of that the unknown parameters are derived from the target eigenpairs, which will be highlighted in Section4. In this paper, the FE analysis (the steps marked with dashed lines in the RMMS) acts as a mostly representative way to generate the target eigenpairs since its correctness and reliability have been recognized, but is not essential if satisfactory target eigeninformation can be obtained from other ways such as experiments.


## 3. Part I: Division and Pre-Evaluation

In contrast with the RFEMs, the proposed RMMS directly divide the flexible structure into RFEs connected by SDEs, which can be combined with mode shapes from FEA to resolve the problem of pre-evaluation of the division scheme. This division step into subdomains can be summarized as a block-wise local-rigid process, which means that the kinematic motion of each subdomain is synchronous. As shown in Figure2, a flexible structure is divided into several synchronistic subdomains or RFEs under a certain division scheme.Figure 2
Division of a flexible structure into RFEs.The generalized displacement vector of the centroid of theith RFE with respect to global frame is defined as(1)qi=riTθiTT=xiyiziαiβiγiT,where ri=xiyiziT are translational displacements and θi=αiβiγiT are the rotational displacements represented by Cardan angles.When assuming that rotational displacements are small enough, displacements of nodej in the ith RFE, denoted as uij, can be expressed as [38](2)uij=Tijqi,where Tij is called the displacement transformation matrix, and the detailed expression is(3)Tij=I−r˜ij0I,where r˜ij is a 3×3 skew symmetric matrix about the vector rij=xijyijzijT of node j relative to the centroid of the ith RFE [38].The total transformation equation of theith RFE can be obtained by assembling all the node-centroid transformations that are expressed as equation (2):(4)Zi=ui1Tui2T⋯uiniTT=Diqi,where ni denotes the total number of nodes in the ith RFE; Zi denotes the displacements of all nodes in the ith RFE; and Di=Ti1TTi2T⋯TiniTT is transformation matrix from the centroid to all nodes in the ith RFE.Assuming that the flexible structure is divided intonb synchronistic subdomains, i.e., RFEs, the generalized displacements of all nodes in the whole flexible structure can be defined as(5)Z=Z1TZ2T⋯ZnbTT=DX,where X=q1Tq2T⋯qnbTT is displacement vector of all centroids of the RFEs and D=diagD1,D2,…,Dnb is total transformation matrix from all centroids to all nodes in the flexible structure.Considering kinematic motions of subdomains under mode shapes in FEA results are asynchronous except for free modes, the MAC [39] is introduced to preliminarily evaluate the synchronization errors between mode shapes in FEA results, and the quasirigid body mode shapes represented by equation (5). Value of MAC presents the similarity of two mode shapes, ranging from 0 to 1. The closer the value is to 1, the higher the similarity between them. The expression of MAC is as follows:(6)MAC=φaTφb2φaTφaφbTφb,where φa and φb are two arbitrary mode shape vectors.DenoteωAk,φAkk=1,2,…,N, the mode frequencies, and mode shapes of the flexible structure in FEA results, and R1,R2,…,Rnb are the nb divided subdomains under a certain division scheme. φRk is the kth mode shape of all the centroids, which fits equation (5), i.e., φAk=DφRk. These linear equations are usually incompatible due to the existence of the synchronization errors. Usually, matrix D is full column rank, and the minimal norm least square solution, utilizing the Moore–Penrose generalized inverse, can be calculated as(7)φRk=D+φAk=DTD−1DTφAk.Equation (7) gives the solution of φRk with a minimum synchronization error under a certain division scheme, and then the quasirigid body mode shape φ⌣k=DφRk=DDTD−1DTφAk. Hence, the MAC value of mode shape vectors φ⌣k and φAk is simplified as(8)MACRk=φ⌣kTφAk2φ⌣kTφ⌣kφAkTφAk=φAkTDDTD−1DTφAk.It is noteworthy that equation (7) gives ideal mode shapes of a rigid MBS under a certain division scheme. However, the final equivalent MBS after the parameter identification usually cannot keep these mode shapes, which means the final MAC value is always smaller than the result obtained from equation (8). For this reason, equation (8) gives the lower bound of the division error range, which can be used to preliminarily evaluate whether the division scheme is reasonable.


## 4. Part II: The Parameter Identification Problem

Once the pre-evaluation is accomplished, results calculated by means of equation (7) can be regarded as target mode shapes of the equivalent rigid MBS, and the corresponding mode frequencies can also be obtained from FEA results. The core problem to be solved is how to determine the dynamic parameters of SDEs, which is actually an inverse eigenvalue problem and can be elaborated as follows: given the mode frequencies ω1<ω2<⋯<ωt and mode shapes φ1,φ2,…,φt, finding the parameters pii=1,2,…,np such that eigenvalues and eigenvectors of the object system are exactly λ1<λ2<⋯<λtλi=ωi2 and φ1,φ2,…,φt, respectively.When considering that both the dynamic parameters and pose parameters are unknown, this problem will be intricate as equations between the system stiffness matrix and variablespi become strongly nonlinear. To avoid this, we firstly prove that the effects of SDEs at different positions and orientations between two RFEs can be equaled. Then, by rationally choosing independent variables, the nonlinear equations can be transformed into linear ones in Rn×n matrix space; in other words, the stiffness matrix can be written as(9)Kp1,p2,…,pnp=K0+∑piKi,where Ki are constant matrices.Furthermore, by using the vectorization operator and Kronecker product [40], linear matrix function in Rn×n space is transformed into a linear combination of vectors in Rn2 space, and the inverse eigenvalue problem is eventually transferred into a least square solution problem of linear equations in Rn2 vector space.


### 4.1. The Equivalence between SDEs

As shown in Figure3, Bi and Bj are two RFEs, Ci and Cj are the centroids, respectively, while qi and qj are the generalized displacements of them. These two RFEs are connected by an SDE at point P (denoted as SDE-P) in initial state, and this section investigates whether it can be equaled to point Q.Figure 3
A schematic diagram of SDEs at different positions and orientations between two RFEs.Detailed definition of the generalized displacement of each centroid can be obtained from equation (1). When α,β,andγ are small angles, the direction cosine matrix of global frame to local frame can be simplified as(10)A≈1−γβγ1−α−βα1=I+θ˜,where θ˜ is also a skew symmetric matrix of size 3×3 about vector θ.DenoteRrb=diagArb,Arb; the generalized displacement with respect to different reference frames can be uniformly expressed as [38](11)qi=Rrbqib.Based on equations (2) and (11), deformation of the SDE-P with respect to its local frame takes the following form:(12)δijPP=RrkPTTiPqi−TjPqj,where RrkP defines the composite direction cosine matrix of SDE-P with respect to the global frame.k
i
j
P denotes the stiffness matrix of SDE-P with respect to local frame {P}, thus the generalized force containing both the force and moment can be written as(13)fijP=FijPtijP=RrkPkijPRrkPTTiPqi−TjPqj.According to the basic law on parallel translation of force [41], force of SDE-P can be shifted to point Q, and the equivalent generalized force is(14)f⌣ijQ=F⌣ijQt⌣ijQ=FijPtijP+rQP×FijP=TQPTfijP,where TQPT=I−r˜QP0IT.Substituting equation (13) into equation (14) yields(15)f⌣ijQ=TQPTRrkPkijPRrkPTTiPqi−TjPqj.With consideration of the transitivity and commutativity of the displacement transformation matrix,TiP=TiQTQP=TQPTiQ, equation (15) can be rewritten as(16)f⌣ijQ=TQPTRrkPkijPRrkPTTQPTiQqi−TjQqj.Assuming that an SDE, whose stiffness matrix iskijQ, exists at point Q (denoted as SDE-Q) between RFE Bi and Bj, the same with equation (13), the generalized force of the SDE-Q is(17)fijQ=FijQtijQ=RrkQkijQRrkQTTiQqi−TjQqj.The relationship between stiffness matriceskijP and kijQ can be obtained from the equation fijQ=f⌣ijQ:(18)kijQ=RrkQTTQPTRrkPkijPRrkPTTQPRrkQ.In this formula, matricesR and T are definitely determined by the relative position vector rPQ and the two reference frames defining stiffness matrices of SDEs and are irrelevant to state variables such as generalized displacements of centroids. Thus, this derivation process proves that an SDE at a certain position and orientation between two RFEs can be equaled to a specified position and orientation, which implies that pose parameters of an SDE are not independent with dynamic parameters of the SDE. In the next section, it can be seen that this proof contributes to the decoupling and linearization.


### 4.2. Decoupling and Linearization

In this section, by choosing appropriate independent variables and introducing the vectorization operator, we will build a linear relationship between the global stiffness matrix and unknown parameters of SDEs inRn2 linear space. As shown in Figure 3, the SDE-Q between RFEs Bi and Bj is taken as an example. The same with equtaion (12), deformation of SDE-Q is(19)δijQQ=RrkQTTiQqi−TjQqj,in which, the definition of RrkQ can refer to equations (11) and (12).The elastic potential energy of SDE-Q is(20)Uij=12δijQQTkijQδijQQ=12qiTTiQT−qjTTjQTRrkQkijQRrkQTTiQqi−TjQqj.The elastic potential energy of the equivalent MBS is(21)Uk=∑i=1nb−1∑j=i+1nbUij,where nb is the total number of the RFEs of the equivalent MBS.The derivative of equation (20) takes the form(22)∂Uij∂qi=TiQTRrkQkijQRrkQTTiQqi−TiQTRrkQkijQRrkQTTjQqj.This formula can be rewritten as(23)∂Uij∂qi=Kiiqi+Kijqj,where(24)Kii=TiQTRrkQkijQRrkQTTiQ,Kij=−TiQTRrkQkijQRrkQTTjQ.Similarly,(25)∂Uij∂qj=Kjiqi+Kjjqj,where(26)Kji=KijT,Kjj=TjQTRrkQkijQRrkQTTjQ.The contribution of SDE between RFEsi and j to the global stiffness matrix is outlined in Figure 4. Extract items related to the SDE between RFEs i and j from the global stiffness matrix and fill the other items with 0 to keep the matrix order, and then denote it as K^ij. Accordingly, the global stiffness matrix can be expressed as(27)KR=∑i=1nb−1∑j=i+1nbK^ij.Figure 4
Contribution of SDE between RFEi and RFE j to the global stiffness matrix.The stiffness matrix of SDE-Q, which is symmetric positive definite, can be partitioned askijQ=kTkTθkTθTkθ. Substituting the direction transformation matrix RrkQ=diagA,A and the displacement transformation matrix TuQ=I−r˜uQ0Iu=i,j, matrix Kuvu,v=i,j takes the form(28)Kuv=δuvAkTATAkTθAT−AkTATr˜vQr˜uQAkTAT+AkθTATAkθAT+r˜uQAkTθAT−AkθTATr˜vQ−r˜uQAkTATr˜vQ,where A=ArkQ, and riQ and rjQ are vectors from centroids of bodies i and j to point Q, respectively, and δuv=1u=v−1u≠v is a sign function.In equation (28), there exist 21 unknown variables describing the stiffness matrix and 6 unknown variables describing the position and orientation of SDE-Q. These 27 unknown variables are not independent with each other, and unreasonable choosing of independent variables will complicate the problem. Nevertheless, referring to the proof in the previous chapter, the decoupling can be accomplished by specifying the position and orientation of SDE-Q. Hence, we choose the 21 unknown parameters in stiffness matrix as independent variables, and the remaining parameters are defined as follows:(1)
Position of SDE-Q is defined as the midpoint of the two centroids; then,riQ=rij/2 and rjQ=−rij/2(2)
Orientation of SDE-Q is defined as parallel to the global frame; then,A=ArkQ=ISincerij is a known vector, linear relationship exists between the global matrix and the independent variables in Rn×n matrix space, namely, the global matrix in equation (27) can be expressed as equation (9).After these decoupling processes, the vectorization operator is utilized to transfer matrices inRn×n matrix space into vectors in Rn2 linear space, advantages of which will be seen in subsequent derivations.DenotekUij as the upper triangular part of kijQ and VeckUij the vectorization [40], the size of which is 21×1. With consideration of the property of Kronecker product [40], vectorization of equations (24) and (26) is(29)VecKii=VecTiQTkTiQ=TiQT⊗TiQTT¯UVeckUij,where T¯U is a constant transformation matrix of size 36×21, between which fits VeckijQ=T¯UVeckUij.Similarly,(30)VecKij=−TjQT⊗TiQTT¯UVeckUij,VecKji=−TiQT⊗TjQTT¯UVeckUij,VecKjj=TjQT⊗TjQTT¯UVeckUij.In the same way, a fixed transformation relation exist between vectorization ofK^ij with VecKuv,u,v=i,j, and this can be denoted as T¯ij, which is directly determined by body numbers i and j. VecK^ij can be written as(31)VecK^ij=T¯ijVecKiiTVecKijTVecKjiTVecKjjTT.Substituting equations (29) and (30) into equation (31) yields(32)VecK^ij=T¯ijTiQT⊗TiQT−TjQT⊗TiQT−TiQT⊗TjQTTjQT⊗TjQTT¯UVeckUij=T¯UijVeckUij,where T¯Uij can be regarded as the final transformation matrix between the global matrix and independent variables of SDE-Q connecting RFEs i and j.In terms of the whole equivalent MBS, combining all the transformation relationships described in equation (32), vectorization of the global matrix takes the form(33)VecKR=∑i=1nb−1∑j=i+1nbVecK^ij=T¯U11⋯T¯Uij⋯VeckU11⋮VeckUij⋮=Wupuk,where puk=VeckU11T⋯VeckUijT⋯T is the final vector of all unknown variables, the size of which is 21nb×1 in SDEs and Wu=T¯U11⋯T¯Uij⋯T is the final transformation matrix between the VecKR and puk, the size of which is 6nb2×21nb when all 6 DOFs of RFEs are considered.


### 4.3. Solving the Parameter Identification Problem

Since damping coefficients of common structures are usually small and have little influence on the modal analysis, theyare ignored in this paper. The equation of motion for the free vibration of the equivalent MBS withn=6nb DOFs constructed by the RMMS can be expressed by(34)MRX¨+KRX=0,where X=q1Tq2T⋯qnbTT is the vector of the generalized displacements defined as in equations (1) and (5), MR is the known global mass matrix, which is symmetric positive definite, and KR is the global stiffness matrix to be solved, which is symmetric positive definite or semidefinite.The associated generalized eigenvalue problem of this MBS is(35)KR−λMRφ=0,where λk=ωk2.The Cholesky decomposition ofMR is MR=LLT, where L is a lower triangular matrix. Substituting this into equation (35) and then premultiplying the equation by L−1 gives(36)L−1KRL−1TLTφ=λLTφ.LetB=L−1KRL−1T and v=LTφ, thus equation (35) is transformed into standard eigenvalue problem Bv=λv. Since L is a known matrix, the corresponding generalized inverse eigenvalue problem can be transformed into a standard inverse eigenvalue problem.λ
=
λ
1
λ
2
⋯
λ
n
T and V=v1v2⋯vn denote the target eigenvalues and eigenvectors, respectively. Thus, we need to solve matrix B under the constraints that(37)Bvi=λivi,i=1,2,…,n.Owing to the orthogonal of eigenvectors, if not, then the eigenvectors should be corrected as in reference [42], and equation (37) can be rewritten as(38)B=VΛVT,where Λ=diagλ1,λ2,…,λn.Applying the vectorization operator to equation (38) gives(39)VecB=V⊗VVecΛ.Substituting equations (33) and (36) into equation (39) yields(40)L−1⊗L−1Wupuk=HVecΛ,where H=V⊗V, Wu and puk are expressed as equation (33).Any column of matrixH can be expressed as VecvsvtTs,t∈1,2,…,n, and by extracting two arbitrary columns Hi=VecvisvitT and Hj=VecvjsvjtT from H, we obtain(41)HiTHj=vis1vitTvis2vitT⋯visnvitTvjs1vjtvjs2vjt⋮vjsnvjt=visTvjsvitTvjt=1,ifi=j,0,ifi≠j..Therefore, all column vectors of matrix exactly form the standard orthonormal basis ofRn2 linear space. Premultiplying equation (40) by HT gives(42)Wppuk=VecΛ,where Wp=HTL−1⊗L−1Wu.The left item of equation (42) forms a subspace in Rn2 linear space determined by the independent variables vector puk, while the right item is the target vector constructed by target eigenvalues. In fact, the eigenpairs of the final MBS constructed by the RMMS cannot perfectly match the target eigenpairs; thus, equation (42) is an essentially incompatible equation. The least squares solution can be obtained using the Moore–Penrose generalized inverse, and the final solution minimizes the Wppuk−VecΛ.If the abovementioned formula is utilized directly, high order modes will dominate. Also, the low order modes will be seriously inconsistent eventually. To avoid this situation, a weight matrixWt, in which diagonal items represent weight coefficients of each order, is introduced to revise equation (42) as(43)W^ppuk=VecWtΛWt,where W^p=Wt⊗WtHTL−1⊗L−1Wu is the modified coefficient matrix.The least squares solution, obtained by Moore–Penrose generalized inverse, ispuk=W^p+VecWtΛWt. Generally, W^p is generally a column-full rank matrix; thus, the final least square solution is(44)puk=W^pTW^p−1W^pTVecWtΛWt.In summary, by introducing a vectorization operator, the parameter identification of the irregular flexible structure is transformed from the initial inverse eigenvalue problem into solving incompatible linear equations. Optimal stiffness parameters of the equivalent rigid MBS can be calculated by equation (44) by means of adjusting the weight coefficients of different orders. The full computing flow of the parameter identification and the RMMS is summarized as in Figure 5.Figure 5
Full computing flow chart of the RMMS.


## 4.1. The Equivalence between SDEs

As shown in Figure3, Bi and Bj are two RFEs, Ci and Cj are the centroids, respectively, while qi and qj are the generalized displacements of them. These two RFEs are connected by an SDE at point P (denoted as SDE-P) in initial state, and this section investigates whether it can be equaled to point Q.Figure 3
A schematic diagram of SDEs at different positions and orientations between two RFEs.Detailed definition of the generalized displacement of each centroid can be obtained from equation (1). When α,β,andγ are small angles, the direction cosine matrix of global frame to local frame can be simplified as(10)A≈1−γβγ1−α−βα1=I+θ˜,where θ˜ is also a skew symmetric matrix of size 3×3 about vector θ.DenoteRrb=diagArb,Arb; the generalized displacement with respect to different reference frames can be uniformly expressed as [38](11)qi=Rrbqib.Based on equations (2) and (11), deformation of the SDE-P with respect to its local frame takes the following form:(12)δijPP=RrkPTTiPqi−TjPqj,where RrkP defines the composite direction cosine matrix of SDE-P with respect to the global frame.k
i
j
P denotes the stiffness matrix of SDE-P with respect to local frame {P}, thus the generalized force containing both the force and moment can be written as(13)fijP=FijPtijP=RrkPkijPRrkPTTiPqi−TjPqj.According to the basic law on parallel translation of force [41], force of SDE-P can be shifted to point Q, and the equivalent generalized force is(14)f⌣ijQ=F⌣ijQt⌣ijQ=FijPtijP+rQP×FijP=TQPTfijP,where TQPT=I−r˜QP0IT.Substituting equation (13) into equation (14) yields(15)f⌣ijQ=TQPTRrkPkijPRrkPTTiPqi−TjPqj.With consideration of the transitivity and commutativity of the displacement transformation matrix,TiP=TiQTQP=TQPTiQ, equation (15) can be rewritten as(16)f⌣ijQ=TQPTRrkPkijPRrkPTTQPTiQqi−TjQqj.Assuming that an SDE, whose stiffness matrix iskijQ, exists at point Q (denoted as SDE-Q) between RFE Bi and Bj, the same with equation (13), the generalized force of the SDE-Q is(17)fijQ=FijQtijQ=RrkQkijQRrkQTTiQqi−TjQqj.The relationship between stiffness matriceskijP and kijQ can be obtained from the equation fijQ=f⌣ijQ:(18)kijQ=RrkQTTQPTRrkPkijPRrkPTTQPRrkQ.In this formula, matricesR and T are definitely determined by the relative position vector rPQ and the two reference frames defining stiffness matrices of SDEs and are irrelevant to state variables such as generalized displacements of centroids. Thus, this derivation process proves that an SDE at a certain position and orientation between two RFEs can be equaled to a specified position and orientation, which implies that pose parameters of an SDE are not independent with dynamic parameters of the SDE. In the next section, it can be seen that this proof contributes to the decoupling and linearization.


## 4.2. Decoupling and Linearization

In this section, by choosing appropriate independent variables and introducing the vectorization operator, we will build a linear relationship between the global stiffness matrix and unknown parameters of SDEs inRn2 linear space. As shown in Figure 3, the SDE-Q between RFEs Bi and Bj is taken as an example. The same with equtaion (12), deformation of SDE-Q is(19)δijQQ=RrkQTTiQqi−TjQqj,in which, the definition of RrkQ can refer to equations (11) and (12).The elastic potential energy of SDE-Q is(20)Uij=12δijQQTkijQδijQQ=12qiTTiQT−qjTTjQTRrkQkijQRrkQTTiQqi−TjQqj.The elastic potential energy of the equivalent MBS is(21)Uk=∑i=1nb−1∑j=i+1nbUij,where nb is the total number of the RFEs of the equivalent MBS.The derivative of equation (20) takes the form(22)∂Uij∂qi=TiQTRrkQkijQRrkQTTiQqi−TiQTRrkQkijQRrkQTTjQqj.This formula can be rewritten as(23)∂Uij∂qi=Kiiqi+Kijqj,where(24)Kii=TiQTRrkQkijQRrkQTTiQ,Kij=−TiQTRrkQkijQRrkQTTjQ.Similarly,(25)∂Uij∂qj=Kjiqi+Kjjqj,where(26)Kji=KijT,Kjj=TjQTRrkQkijQRrkQTTjQ.The contribution of SDE between RFEsi and j to the global stiffness matrix is outlined in Figure 4. Extract items related to the SDE between RFEs i and j from the global stiffness matrix and fill the other items with 0 to keep the matrix order, and then denote it as K^ij. Accordingly, the global stiffness matrix can be expressed as(27)KR=∑i=1nb−1∑j=i+1nbK^ij.Figure 4
Contribution of SDE between RFEi and RFE j to the global stiffness matrix.The stiffness matrix of SDE-Q, which is symmetric positive definite, can be partitioned askijQ=kTkTθkTθTkθ. Substituting the direction transformation matrix RrkQ=diagA,A and the displacement transformation matrix TuQ=I−r˜uQ0Iu=i,j, matrix Kuvu,v=i,j takes the form(28)Kuv=δuvAkTATAkTθAT−AkTATr˜vQr˜uQAkTAT+AkθTATAkθAT+r˜uQAkTθAT−AkθTATr˜vQ−r˜uQAkTATr˜vQ,where A=ArkQ, and riQ and rjQ are vectors from centroids of bodies i and j to point Q, respectively, and δuv=1u=v−1u≠v is a sign function.In equation (28), there exist 21 unknown variables describing the stiffness matrix and 6 unknown variables describing the position and orientation of SDE-Q. These 27 unknown variables are not independent with each other, and unreasonable choosing of independent variables will complicate the problem. Nevertheless, referring to the proof in the previous chapter, the decoupling can be accomplished by specifying the position and orientation of SDE-Q. Hence, we choose the 21 unknown parameters in stiffness matrix as independent variables, and the remaining parameters are defined as follows:(1)
Position of SDE-Q is defined as the midpoint of the two centroids; then,riQ=rij/2 and rjQ=−rij/2(2)
Orientation of SDE-Q is defined as parallel to the global frame; then,A=ArkQ=ISincerij is a known vector, linear relationship exists between the global matrix and the independent variables in Rn×n matrix space, namely, the global matrix in equation (27) can be expressed as equation (9).After these decoupling processes, the vectorization operator is utilized to transfer matrices inRn×n matrix space into vectors in Rn2 linear space, advantages of which will be seen in subsequent derivations.DenotekUij as the upper triangular part of kijQ and VeckUij the vectorization [40], the size of which is 21×1. With consideration of the property of Kronecker product [40], vectorization of equations (24) and (26) is(29)VecKii=VecTiQTkTiQ=TiQT⊗TiQTT¯UVeckUij,where T¯U is a constant transformation matrix of size 36×21, between which fits VeckijQ=T¯UVeckUij.Similarly,(30)VecKij=−TjQT⊗TiQTT¯UVeckUij,VecKji=−TiQT⊗TjQTT¯UVeckUij,VecKjj=TjQT⊗TjQTT¯UVeckUij.In the same way, a fixed transformation relation exist between vectorization ofK^ij with VecKuv,u,v=i,j, and this can be denoted as T¯ij, which is directly determined by body numbers i and j. VecK^ij can be written as(31)VecK^ij=T¯ijVecKiiTVecKijTVecKjiTVecKjjTT.Substituting equations (29) and (30) into equation (31) yields(32)VecK^ij=T¯ijTiQT⊗TiQT−TjQT⊗TiQT−TiQT⊗TjQTTjQT⊗TjQTT¯UVeckUij=T¯UijVeckUij,where T¯Uij can be regarded as the final transformation matrix between the global matrix and independent variables of SDE-Q connecting RFEs i and j.In terms of the whole equivalent MBS, combining all the transformation relationships described in equation (32), vectorization of the global matrix takes the form(33)VecKR=∑i=1nb−1∑j=i+1nbVecK^ij=T¯U11⋯T¯Uij⋯VeckU11⋮VeckUij⋮=Wupuk,where puk=VeckU11T⋯VeckUijT⋯T is the final vector of all unknown variables, the size of which is 21nb×1 in SDEs and Wu=T¯U11⋯T¯Uij⋯T is the final transformation matrix between the VecKR and puk, the size of which is 6nb2×21nb when all 6 DOFs of RFEs are considered.


## 4.3. Solving the Parameter Identification Problem

Since damping coefficients of common structures are usually small and have little influence on the modal analysis, theyare ignored in this paper. The equation of motion for the free vibration of the equivalent MBS withn=6nb DOFs constructed by the RMMS can be expressed by(34)MRX¨+KRX=0,where X=q1Tq2T⋯qnbTT is the vector of the generalized displacements defined as in equations (1) and (5), MR is the known global mass matrix, which is symmetric positive definite, and KR is the global stiffness matrix to be solved, which is symmetric positive definite or semidefinite.The associated generalized eigenvalue problem of this MBS is(35)KR−λMRφ=0,where λk=ωk2.The Cholesky decomposition ofMR is MR=LLT, where L is a lower triangular matrix. Substituting this into equation (35) and then premultiplying the equation by L−1 gives(36)L−1KRL−1TLTφ=λLTφ.LetB=L−1KRL−1T and v=LTφ, thus equation (35) is transformed into standard eigenvalue problem Bv=λv. Since L is a known matrix, the corresponding generalized inverse eigenvalue problem can be transformed into a standard inverse eigenvalue problem.λ
=
λ
1
λ
2
⋯
λ
n
T and V=v1v2⋯vn denote the target eigenvalues and eigenvectors, respectively. Thus, we need to solve matrix B under the constraints that(37)Bvi=λivi,i=1,2,…,n.Owing to the orthogonal of eigenvectors, if not, then the eigenvectors should be corrected as in reference [42], and equation (37) can be rewritten as(38)B=VΛVT,where Λ=diagλ1,λ2,…,λn.Applying the vectorization operator to equation (38) gives(39)VecB=V⊗VVecΛ.Substituting equations (33) and (36) into equation (39) yields(40)L−1⊗L−1Wupuk=HVecΛ,where H=V⊗V, Wu and puk are expressed as equation (33).Any column of matrixH can be expressed as VecvsvtTs,t∈1,2,…,n, and by extracting two arbitrary columns Hi=VecvisvitT and Hj=VecvjsvjtT from H, we obtain(41)HiTHj=vis1vitTvis2vitT⋯visnvitTvjs1vjtvjs2vjt⋮vjsnvjt=visTvjsvitTvjt=1,ifi=j,0,ifi≠j..Therefore, all column vectors of matrix exactly form the standard orthonormal basis ofRn2 linear space. Premultiplying equation (40) by HT gives(42)Wppuk=VecΛ,where Wp=HTL−1⊗L−1Wu.The left item of equation (42) forms a subspace in Rn2 linear space determined by the independent variables vector puk, while the right item is the target vector constructed by target eigenvalues. In fact, the eigenpairs of the final MBS constructed by the RMMS cannot perfectly match the target eigenpairs; thus, equation (42) is an essentially incompatible equation. The least squares solution can be obtained using the Moore–Penrose generalized inverse, and the final solution minimizes the Wppuk−VecΛ.If the abovementioned formula is utilized directly, high order modes will dominate. Also, the low order modes will be seriously inconsistent eventually. To avoid this situation, a weight matrixWt, in which diagonal items represent weight coefficients of each order, is introduced to revise equation (42) as(43)W^ppuk=VecWtΛWt,where W^p=Wt⊗WtHTL−1⊗L−1Wu is the modified coefficient matrix.The least squares solution, obtained by Moore–Penrose generalized inverse, ispuk=W^p+VecWtΛWt. Generally, W^p is generally a column-full rank matrix; thus, the final least square solution is(44)puk=W^pTW^p−1W^pTVecWtΛWt.In summary, by introducing a vectorization operator, the parameter identification of the irregular flexible structure is transformed from the initial inverse eigenvalue problem into solving incompatible linear equations. Optimal stiffness parameters of the equivalent rigid MBS can be calculated by equation (44) by means of adjusting the weight coefficients of different orders. The full computing flow of the parameter identification and the RMMS is summarized as in Figure 5.Figure 5
Full computing flow chart of the RMMS.


## 5. Experimental Verification and Application

In this section, two cases were carried out. Among them, Case I was mainly designed to validate the feasibility and correctness of the parameter identification method and the RMMS; experiments were conducted on a planar motor mover with Halbach magnet array, and the simulated modal behaviors from models established by the RMMS were compared with the FEA ones and the experimental results, respectively. Case II was mainly conducted to exhibit the applicability of the RMMS for large-scale mechanical systems; the RMMS was applied in dynamic modeling of a lithography system.


### 5.1. Case I: Experimental Verification on a Planar Motor

#### 5.1.1. Set-Up of Experiments

As planar motor widely applied in high-precision equipment, the inherent structural characteristics of the mover must be considered to meet the requirements for near-nano precision control. The employed planar motor mover, as depicted in Figure6, mainly comprises two parts: an aluminum frame and a magnets array. The frame has Young’s modulus E=71GPa, material density ρ=2.77g/cm3, and Poisson’s ratio μ=0.33. The magnets array is made up of horizontal and vertical magnetized N45 permanent magnets, which have E=160GPa, ρ=7.4g/cm3, μ=0.24. The fundamental structural parameters of the planar motor mover are listed in Table 1.Figure 6
A schematic of the planar motor mover with Halbach magnet array.Table 1
Structural parameters of the planar motor mover.
Items
Symbol
Value (mm)
Pole pitch
τ
30
Width of magnets
τ
m
20.5
Height of magnets
H
m
7
Width of frame
W
f
360
Height of frame
H
f
15Modal experiments were conducted on a test rig, as shown in Figure7. The planar motor mover was supported by a sponge to ensure that the normal modal results were measured. The external excitation was imposed by an impact hammer (PCB 086C01) with a hard trigger head, and a roving hammer strategy was employed to obtain vibrational behaviors of different nodes. Three PCB 356A24 acceleration sensors were placed at different positions to measure the vibration response. The data acquisition device was a LMS SCADAS Ш, and its accompanying software LMS Test.Lab was utilized for data processing and analysis.Figure 7
Photograph of the test rig.


#### 5.1.2. Description of Simulation Models

In the planar mover, the magnets were assembled by means of bonding with special glue; thus, thin and soft adhesive layers were formed between magnets, as well as between the array and the frame, as shown in Figure8. Due to the heterogeneity and discontinuity of the adhesive layers introduced by assembly errors and other factors, it is extraordinarily difficult to establish a FE model which can accurately predict the vibrational behaviors, and at least, the resulted model order is very high.Figure 8
Partial section views of the adhesive layer in the planar motor mover.To illustrate the feasibility and advantages of the proposed RMMS from different perspectives, this paper constructed three kinds of simulation models based on the proposed method and the FEA method. The FE model was performed with the ANSYS Workbench, while the RMMS models with MATLAB platform. They are described and denoted as follows:(1)
FEM: a rough FE model which ignores the effects of the adhesive layers, the magnet array is modeled as one part with N45 materials, and nodes at the interface between it and the aluminum frame are boned(2)
MBSFm_n (without modification): an equivalent MBS constructed by the proposed method, which divides the planar mover intom×n rigid bodies, and the parameters are derived from the FEA results directly(3)
MBSMm_n (with modification): an equivalent MBS constructed by the proposed method, which divides the planar mover intom×n rigid bodies, and the parameters are derived from the FEA results that are modified by the experimental results.Figure9 shows the schematic of the simulation models, including the established MBS with division scheme (m×n) and the FE model. In the following, firstly, both the static and the modal behaviors analyzed by MBSFs with no modification in parameter identification are compared with FEA results. Secondly, modal results by MBSMs with experimental modification in parameter identification are compared with the experimental ones. As mentioned before, the damping coefficients of common structures are usually small and ignored during the parameter identification in RMMS. For comparison purposes on modal behaviors, the modal damping ratios of each mode in simulation results are set to be consistent with that of the relevant experimental modes, which are about 5%.Figure 9
Schematic of the simulation models: (a) equivalent MBS with division scheme (m×n) and (b) FE model.
(a)
(b)


#### 5.1.3. Comparison of MBSFs Results and FEA Results

In the static analysis, outer region of the aluminum frame (marked with light green in Figure9(b)) was fixed and a 100 N force was applied on the central point (area). Correspondingly, the outer RFEs in MBSFs were fixed, and a same force with 100 N was loaded on the central RFE. The force was loaded step by step, and deformations of the central point computed by the two type models are compared in Figure 10. It can be seen that the force-displacement curves are linear in the analysis due to the high structural stiffness. Meanwhile, the central deformations computed by the MBSFs are in good accordance with the FEA results as the division density is highly enough.Figure 10
Static analysis results of the planar motor mover: (a) FEA results and (b) the MBSFs results.
(a)
(b)Modal behaviors were analyzed by both the constructed MBSFs and FE model on a desktop computer that processes a 3.2 GHz Intel Core i5 CPU and 16 GB of RAM. Both the model scales and the computing times are compared in Table2. The results show that the RMMS can dramatically reduce the model scale since model DOFs of the two MBSFs are much lower than that of the FE model, and thus greatly shortens the computing time.Table 2
Model scales and computing times of the MBSFs and FE model.
Items
MBSFs
FEM
MBSF3_3
MBSF10_10
Rigid bodies/nodes
9
100
53381
Model DOFs
54
600
160143
Computing time (s)
0.00051
0.042
18.41Figure11 presents the MAC matrices that show mode shape similarity of MBSF3_3 vs. FEA modes and MBSF10_10 vs. FEA modes, respectively. In contrast with that in Figure 11(a), the MAC matrix in Figure 11(b) has high diagonal values that are close to 1, excepting for the 4th and 5th modes, which indicates that the mode shapes correlated well between the two mode sets from the MBSF10_10 and the FE model. It is worth noting that mode 4 and 5 computed by the FE model actually have a repeated eigenvalue for the structural symmetry; thus, the related eigenvectors are not unique and hard to be matched.Figure 11
MAC matrix showing mode shape similarity between modes from MBSF models and the FE model: (a) MBSF3_3 vs. FEA and (b) MBSF10_10 vs. FEA.
(a)
(b)Figure11 also implies the correspondence between the MBSFs and FEA modes, which connects the similar mode-pairs together. Based on this correlation, frequency relative errors of the first 10 nonzero modes with respect to the FEA results are calculated and listed in Table 3. As can be seen that the maximum absolute value of the relative errors between the MBSF10_10 and the FEA modes is less than 5%, the average value is 3.4%.Table 3
Relative errors of mode frequencies with respect to FEA results.
Mode
Frequency (Hz)
Absolute value of relative error (%)
FEA
MBSF3_3
MBSF10_10
MBSF3_3
MBSF10_10
1
377.4
388.7
394.7
3.0
4.6
2
556.1
568.6
588.0
2.2
5.7
3
674.2
618.0
644.9
8.3
4.4
4
964.8
823.1
976.5
14.7
1.2
5
964.9
943.3
990.5
2.2
2.7
6
1662.6
1250.3
1628.0
24.8
2.1
7
1662.6
1310.5
1629.6
21.2
2.0
8
1735.5
1303.0
1717.4
24.9
1.0
9
1911.3
1735.8
2001.0
9.2
4.7
10
2049.9
1520.4
1933.7
25.8
5.7To compare the vibration behaviors of the two types of models comprehensively, the force-acceleration transfer curves fromz-direction of point P1 to z-direction of point P25 are plotted in Figure 12. It can be observed that the transfer function curve calculated by the MBSF10_10 (marked with solid magenta line) is in good agreement with the curve calculated by the FE model (marked with solid blue line). By comparison, a conclusion can be drawn that the vibration behavior of MBSF constructed by the RMMS, while divided reasonably, is consistent with the FE model.Figure 12
Comparison of force-acceleration transfer functions between 200 and 3000 Hz: P1-z to P25-z.


#### 5.1.4. Comparison of MBSMs Results and Experiment Results

This section further compares the experimental modal results with the simulated ones computed by MBSMs, whose parameters are identified from experimentally modified FEA results. The model information of the MBSMs is listed in Table4.Table 4
Model scales of the MBSMs.
Items
MBSM3_3
MBSM6_6
MBSM8_8
MBSM10_10
Rigid bodies
9
36
64
100
Model DOFs
54
216
384
600Figure13 presents the stabilization diagram of the experimental normal modes in the LMS test. The mode frequencies and related shapes can be obtained in a laboratory. The diagram shows that the measured (MEA) results have good consistency and high reliability for comparison.Figure 13
Stabilization diagram of the experimental results.Figure14 presents the MAC matrices showing the mode shape correlation between the modes of the experiments and that of MBSMs with different division schemes, ranging from 3×3 to 10×10. It can be seen that mode shapes of the MBSM with an insufficient dividing scheme, such as MBSM3_3, are poorly correlated to the experimental ones because its MAC matrix has many high off-diagonal values, whereas good consistencies are achieved under reasonable division since MBSM6_6.Figure 14
MAC matrix showing mode shape similarity between modes from MBSM models and experimental results: (a) MBSM3_3 vs. MEA, (b) MBSM6_6 vs. MEA, (c) MBSM8_8 vs. MEA, and (d) MBSM10_10 vs. MEA.
(a)
(b)
(c)
(d)According to the correspondence between the MBSMs and FEA modes, as implied in MAC graphs of Figure14, their frequency relative errors in the first 10 nonzero modes are calculated and plotted in Figure 15. It can be seen that the maximum relative errors of the MBSMs with divisions from 6×6 to 10×10, with respect to the MEA results, are all less than 10%, and most of them are small than 5%.Figure 15
Mode frequency relative errors with respect to experimental results.In order to intuitively compare the vibration behaviors obtained by the simulation models and the experimental data, the force-acceleration transfer curves from different input channels to different out channels are plotted in Figures16 and 17. The results show that the transfer curves computed by the MBSMs are in good agreement with the experimental ones. It can be concluded that the vibration behaviors of MBSFs, which are established by the RMMS with experimental modification in parameter identification, are consistent with the measured results. Besides, the higher the division density, the better the consistency. However, the increase of rigid elements obviously amplifies the computational burden.Figure 16
Comparison of force-acceleration transfer functions between 200 and 1500 Hz: P1-z to P25-z.Figure 17
Comparison of force-acceleration transfer functions between 200 and 1500 Hz: P9-z to P21-z.


### 5.2. Case II: Application in a Lithography System

An additional application in concern with the dynamic modeling of lithography, which is key to the system design, was presented. The employed lithography, as shown in Figure18(a), mainly consists of an illuminator system, a reticle system, a wafer system, a projection system, precision measurement systems, vibration isolation systems, and frames. Both the reticle system and the wafer system are composed of a coarse stage and a fine stage. Among them, the coarse stage is to realize the micron-level high-acceleration motion, while the fine stage is to achieve the nanolevel ultraprecision motion. The illuminator system, reticle system, wafer system, projection system, and precision measurement system are all installed on the main frame. Also, three vibration isolators are used to connect the main frame with the base frame (i.e., the external world), so as to guarantee an ultrastable environment. In addition, there exist many flexible connections to further isolate the vibration, such as between the main frame and the projection lens, between the main frame and the supporting structures of the two motion stages.Figure 18
Structure of the lithography system: (a) schematic diagram and (b) rigid multibody model.
(a)
(b)In modeling of key structures/subsystems such as the illuminator system, the reticle, and wafer stages, the proposed parameter-identification method and RMMS were applied to establish the equivalent MBSs with parameters identified from modal results of related FE models, as shown in Figure19. Finally, a rigid multibody model of the lithography, which has 75 bodies and 127 DOFs, was constructed and presented in Figure 18(b). With this established model, the modal behavior was analyzed and compared with the experimental results.Figure 19
Dynamic models of key components/subsystems: (a) FEA models and (b) equivalent MBSs.The simulation and experimental modal results within 200 Hz are compared in Table5; it can be observed that the maximum relative error of mode frequency is about 6% while the mode shapes are matched. In addition, Figure 20 compares the simulation and experimental mode frequencies within 500 Hz. It can be summarized that vibration behavior computed by the final established rigid multibody model shows good agreement with the results from measured data, and only small errors exist between measured and the corresponding computational mode orders, especially within 0∼200 Hz, which is a mainly concerned frequency band in the design of lithography. It should also be noted that some modes are not caught in experiments due to the arrangement of sensors and limits of excitations.Table 5
Comparison of simulation and experimental modal results within 200 Hz.
Mode
Mode freq. (Hz)
Error (%)
Mode shape
Simulation
Experiment
Object
Mode description
1
1.49
1.5
0.4
Main frame
Translation alongY
2
1.5
1.51
0.5
Main frame
Translation alongX
3
1.99
2
0.1
Main frame
Translation alongZ
4
2.59
2.61
0.8
Main frame
Rotation aroundX
5
2.67
2.62
1.9
Main frame
Rotation aroundY
6
2.98
3
0.1
Main frame
Rotation aroundZ
7
37.06
39.01
5
Support plate of reticle stage
Torsional vibration
8
39.39
39.03
1
Support plate of reticle stage
Torsional vibration
9
39.47
39.0
1.3
Support plate of reticle stage
Torsional vibration
10
49.7
50
0.6
Wafer stage
Translation alongX
11
69.19
70
1.2
Wafer stage
Translation alongX + rotation around Z
12
76.87
78
1.5
Wafer stage
Rotation aroundX
13
80.33
80
0.5
Reticle stage
Translation alongY
14
82.52
81
1.9
Wafer stage
Rotation aroundZ + translation along X
15
101.9
98
4
Reticle stage
Translation alongX
16
112.9
120
6
Wafer stage
Rotation aroundZ
17
124
128
3.2
Reticle stage
Rotation aroundX
18
127.3
134
5
Lens + reticle stage
—
19
128.4
—
—
Reticle stage
Translation alongZ
20
147.8
150
1.4
Reticle stage
Rotation aroundX
21
171.2
173
1
Lithography
—
22
176.8
182
2.7
Reticle stage
Rotation aroundZ
23
190.8
195
2.2
Reticle stage
Rotation aroundYFigure 20
Comparison of simulation and experimental modal frequencies within 500 Hz.


## 5.1. Case I: Experimental Verification on a Planar Motor

### 5.1.1. Set-Up of Experiments

As planar motor widely applied in high-precision equipment, the inherent structural characteristics of the mover must be considered to meet the requirements for near-nano precision control. The employed planar motor mover, as depicted in Figure6, mainly comprises two parts: an aluminum frame and a magnets array. The frame has Young’s modulus E=71GPa, material density ρ=2.77g/cm3, and Poisson’s ratio μ=0.33. The magnets array is made up of horizontal and vertical magnetized N45 permanent magnets, which have E=160GPa, ρ=7.4g/cm3, μ=0.24. The fundamental structural parameters of the planar motor mover are listed in Table 1.Figure 6
A schematic of the planar motor mover with Halbach magnet array.Table 1
Structural parameters of the planar motor mover.
Items
Symbol
Value (mm)
Pole pitch
τ
30
Width of magnets
τ
m
20.5
Height of magnets
H
m
7
Width of frame
W
f
360
Height of frame
H
f
15Modal experiments were conducted on a test rig, as shown in Figure7. The planar motor mover was supported by a sponge to ensure that the normal modal results were measured. The external excitation was imposed by an impact hammer (PCB 086C01) with a hard trigger head, and a roving hammer strategy was employed to obtain vibrational behaviors of different nodes. Three PCB 356A24 acceleration sensors were placed at different positions to measure the vibration response. The data acquisition device was a LMS SCADAS Ш, and its accompanying software LMS Test.Lab was utilized for data processing and analysis.Figure 7
Photograph of the test rig.


### 5.1.2. Description of Simulation Models

In the planar mover, the magnets were assembled by means of bonding with special glue; thus, thin and soft adhesive layers were formed between magnets, as well as between the array and the frame, as shown in Figure8. Due to the heterogeneity and discontinuity of the adhesive layers introduced by assembly errors and other factors, it is extraordinarily difficult to establish a FE model which can accurately predict the vibrational behaviors, and at least, the resulted model order is very high.Figure 8
Partial section views of the adhesive layer in the planar motor mover.To illustrate the feasibility and advantages of the proposed RMMS from different perspectives, this paper constructed three kinds of simulation models based on the proposed method and the FEA method. The FE model was performed with the ANSYS Workbench, while the RMMS models with MATLAB platform. They are described and denoted as follows:(1)
FEM: a rough FE model which ignores the effects of the adhesive layers, the magnet array is modeled as one part with N45 materials, and nodes at the interface between it and the aluminum frame are boned(2)
MBSFm_n (without modification): an equivalent MBS constructed by the proposed method, which divides the planar mover intom×n rigid bodies, and the parameters are derived from the FEA results directly(3)
MBSMm_n (with modification): an equivalent MBS constructed by the proposed method, which divides the planar mover intom×n rigid bodies, and the parameters are derived from the FEA results that are modified by the experimental results.Figure9 shows the schematic of the simulation models, including the established MBS with division scheme (m×n) and the FE model. In the following, firstly, both the static and the modal behaviors analyzed by MBSFs with no modification in parameter identification are compared with FEA results. Secondly, modal results by MBSMs with experimental modification in parameter identification are compared with the experimental ones. As mentioned before, the damping coefficients of common structures are usually small and ignored during the parameter identification in RMMS. For comparison purposes on modal behaviors, the modal damping ratios of each mode in simulation results are set to be consistent with that of the relevant experimental modes, which are about 5%.Figure 9
Schematic of the simulation models: (a) equivalent MBS with division scheme (m×n) and (b) FE model.
(a)
(b)


### 5.1.3. Comparison of MBSFs Results and FEA Results

In the static analysis, outer region of the aluminum frame (marked with light green in Figure9(b)) was fixed and a 100 N force was applied on the central point (area). Correspondingly, the outer RFEs in MBSFs were fixed, and a same force with 100 N was loaded on the central RFE. The force was loaded step by step, and deformations of the central point computed by the two type models are compared in Figure 10. It can be seen that the force-displacement curves are linear in the analysis due to the high structural stiffness. Meanwhile, the central deformations computed by the MBSFs are in good accordance with the FEA results as the division density is highly enough.Figure 10
Static analysis results of the planar motor mover: (a) FEA results and (b) the MBSFs results.
(a)
(b)Modal behaviors were analyzed by both the constructed MBSFs and FE model on a desktop computer that processes a 3.2 GHz Intel Core i5 CPU and 16 GB of RAM. Both the model scales and the computing times are compared in Table2. The results show that the RMMS can dramatically reduce the model scale since model DOFs of the two MBSFs are much lower than that of the FE model, and thus greatly shortens the computing time.Table 2
Model scales and computing times of the MBSFs and FE model.
Items
MBSFs
FEM
MBSF3_3
MBSF10_10
Rigid bodies/nodes
9
100
53381
Model DOFs
54
600
160143
Computing time (s)
0.00051
0.042
18.41Figure11 presents the MAC matrices that show mode shape similarity of MBSF3_3 vs. FEA modes and MBSF10_10 vs. FEA modes, respectively. In contrast with that in Figure 11(a), the MAC matrix in Figure 11(b) has high diagonal values that are close to 1, excepting for the 4th and 5th modes, which indicates that the mode shapes correlated well between the two mode sets from the MBSF10_10 and the FE model. It is worth noting that mode 4 and 5 computed by the FE model actually have a repeated eigenvalue for the structural symmetry; thus, the related eigenvectors are not unique and hard to be matched.Figure 11
MAC matrix showing mode shape similarity between modes from MBSF models and the FE model: (a) MBSF3_3 vs. FEA and (b) MBSF10_10 vs. FEA.
(a)
(b)Figure11 also implies the correspondence between the MBSFs and FEA modes, which connects the similar mode-pairs together. Based on this correlation, frequency relative errors of the first 10 nonzero modes with respect to the FEA results are calculated and listed in Table 3. As can be seen that the maximum absolute value of the relative errors between the MBSF10_10 and the FEA modes is less than 5%, the average value is 3.4%.Table 3
Relative errors of mode frequencies with respect to FEA results.
Mode
Frequency (Hz)
Absolute value of relative error (%)
FEA
MBSF3_3
MBSF10_10
MBSF3_3
MBSF10_10
1
377.4
388.7
394.7
3.0
4.6
2
556.1
568.6
588.0
2.2
5.7
3
674.2
618.0
644.9
8.3
4.4
4
964.8
823.1
976.5
14.7
1.2
5
964.9
943.3
990.5
2.2
2.7
6
1662.6
1250.3
1628.0
24.8
2.1
7
1662.6
1310.5
1629.6
21.2
2.0
8
1735.5
1303.0
1717.4
24.9
1.0
9
1911.3
1735.8
2001.0
9.2
4.7
10
2049.9
1520.4
1933.7
25.8
5.7To compare the vibration behaviors of the two types of models comprehensively, the force-acceleration transfer curves fromz-direction of point P1 to z-direction of point P25 are plotted in Figure 12. It can be observed that the transfer function curve calculated by the MBSF10_10 (marked with solid magenta line) is in good agreement with the curve calculated by the FE model (marked with solid blue line). By comparison, a conclusion can be drawn that the vibration behavior of MBSF constructed by the RMMS, while divided reasonably, is consistent with the FE model.Figure 12
Comparison of force-acceleration transfer functions between 200 and 3000 Hz: P1-z to P25-z.


### 5.1.4. Comparison of MBSMs Results and Experiment Results

This section further compares the experimental modal results with the simulated ones computed by MBSMs, whose parameters are identified from experimentally modified FEA results. The model information of the MBSMs is listed in Table4.Table 4
Model scales of the MBSMs.
Items
MBSM3_3
MBSM6_6
MBSM8_8
MBSM10_10
Rigid bodies
9
36
64
100
Model DOFs
54
216
384
600Figure13 presents the stabilization diagram of the experimental normal modes in the LMS test. The mode frequencies and related shapes can be obtained in a laboratory. The diagram shows that the measured (MEA) results have good consistency and high reliability for comparison.Figure 13
Stabilization diagram of the experimental results.Figure14 presents the MAC matrices showing the mode shape correlation between the modes of the experiments and that of MBSMs with different division schemes, ranging from 3×3 to 10×10. It can be seen that mode shapes of the MBSM with an insufficient dividing scheme, such as MBSM3_3, are poorly correlated to the experimental ones because its MAC matrix has many high off-diagonal values, whereas good consistencies are achieved under reasonable division since MBSM6_6.Figure 14
MAC matrix showing mode shape similarity between modes from MBSM models and experimental results: (a) MBSM3_3 vs. MEA, (b) MBSM6_6 vs. MEA, (c) MBSM8_8 vs. MEA, and (d) MBSM10_10 vs. MEA.
(a)
(b)
(c)
(d)According to the correspondence between the MBSMs and FEA modes, as implied in MAC graphs of Figure14, their frequency relative errors in the first 10 nonzero modes are calculated and plotted in Figure 15. It can be seen that the maximum relative errors of the MBSMs with divisions from 6×6 to 10×10, with respect to the MEA results, are all less than 10%, and most of them are small than 5%.Figure 15
Mode frequency relative errors with respect to experimental results.In order to intuitively compare the vibration behaviors obtained by the simulation models and the experimental data, the force-acceleration transfer curves from different input channels to different out channels are plotted in Figures16 and 17. The results show that the transfer curves computed by the MBSMs are in good agreement with the experimental ones. It can be concluded that the vibration behaviors of MBSFs, which are established by the RMMS with experimental modification in parameter identification, are consistent with the measured results. Besides, the higher the division density, the better the consistency. However, the increase of rigid elements obviously amplifies the computational burden.Figure 16
Comparison of force-acceleration transfer functions between 200 and 1500 Hz: P1-z to P25-z.Figure 17
Comparison of force-acceleration transfer functions between 200 and 1500 Hz: P9-z to P21-z.


## 5.1.1. Set-Up of Experiments

As planar motor widely applied in high-precision equipment, the inherent structural characteristics of the mover must be considered to meet the requirements for near-nano precision control. The employed planar motor mover, as depicted in Figure6, mainly comprises two parts: an aluminum frame and a magnets array. The frame has Young’s modulus E=71GPa, material density ρ=2.77g/cm3, and Poisson’s ratio μ=0.33. The magnets array is made up of horizontal and vertical magnetized N45 permanent magnets, which have E=160GPa, ρ=7.4g/cm3, μ=0.24. The fundamental structural parameters of the planar motor mover are listed in Table 1.Figure 6
A schematic of the planar motor mover with Halbach magnet array.Table 1
Structural parameters of the planar motor mover.
Items
Symbol
Value (mm)
Pole pitch
τ
30
Width of magnets
τ
m
20.5
Height of magnets
H
m
7
Width of frame
W
f
360
Height of frame
H
f
15Modal experiments were conducted on a test rig, as shown in Figure7. The planar motor mover was supported by a sponge to ensure that the normal modal results were measured. The external excitation was imposed by an impact hammer (PCB 086C01) with a hard trigger head, and a roving hammer strategy was employed to obtain vibrational behaviors of different nodes. Three PCB 356A24 acceleration sensors were placed at different positions to measure the vibration response. The data acquisition device was a LMS SCADAS Ш, and its accompanying software LMS Test.Lab was utilized for data processing and analysis.Figure 7
Photograph of the test rig.


## 5.1.2. Description of Simulation Models

In the planar mover, the magnets were assembled by means of bonding with special glue; thus, thin and soft adhesive layers were formed between magnets, as well as between the array and the frame, as shown in Figure8. Due to the heterogeneity and discontinuity of the adhesive layers introduced by assembly errors and other factors, it is extraordinarily difficult to establish a FE model which can accurately predict the vibrational behaviors, and at least, the resulted model order is very high.Figure 8
Partial section views of the adhesive layer in the planar motor mover.To illustrate the feasibility and advantages of the proposed RMMS from different perspectives, this paper constructed three kinds of simulation models based on the proposed method and the FEA method. The FE model was performed with the ANSYS Workbench, while the RMMS models with MATLAB platform. They are described and denoted as follows:(1)
FEM: a rough FE model which ignores the effects of the adhesive layers, the magnet array is modeled as one part with N45 materials, and nodes at the interface between it and the aluminum frame are boned(2)
MBSFm_n (without modification): an equivalent MBS constructed by the proposed method, which divides the planar mover intom×n rigid bodies, and the parameters are derived from the FEA results directly(3)
MBSMm_n (with modification): an equivalent MBS constructed by the proposed method, which divides the planar mover intom×n rigid bodies, and the parameters are derived from the FEA results that are modified by the experimental results.Figure9 shows the schematic of the simulation models, including the established MBS with division scheme (m×n) and the FE model. In the following, firstly, both the static and the modal behaviors analyzed by MBSFs with no modification in parameter identification are compared with FEA results. Secondly, modal results by MBSMs with experimental modification in parameter identification are compared with the experimental ones. As mentioned before, the damping coefficients of common structures are usually small and ignored during the parameter identification in RMMS. For comparison purposes on modal behaviors, the modal damping ratios of each mode in simulation results are set to be consistent with that of the relevant experimental modes, which are about 5%.Figure 9
Schematic of the simulation models: (a) equivalent MBS with division scheme (m×n) and (b) FE model.
(a)
(b)


## 5.1.3. Comparison of MBSFs Results and FEA Results

In the static analysis, outer region of the aluminum frame (marked with light green in Figure9(b)) was fixed and a 100 N force was applied on the central point (area). Correspondingly, the outer RFEs in MBSFs were fixed, and a same force with 100 N was loaded on the central RFE. The force was loaded step by step, and deformations of the central point computed by the two type models are compared in Figure 10. It can be seen that the force-displacement curves are linear in the analysis due to the high structural stiffness. Meanwhile, the central deformations computed by the MBSFs are in good accordance with the FEA results as the division density is highly enough.Figure 10
Static analysis results of the planar motor mover: (a) FEA results and (b) the MBSFs results.
(a)
(b)Modal behaviors were analyzed by both the constructed MBSFs and FE model on a desktop computer that processes a 3.2 GHz Intel Core i5 CPU and 16 GB of RAM. Both the model scales and the computing times are compared in Table2. The results show that the RMMS can dramatically reduce the model scale since model DOFs of the two MBSFs are much lower than that of the FE model, and thus greatly shortens the computing time.Table 2
Model scales and computing times of the MBSFs and FE model.
Items
MBSFs
FEM
MBSF3_3
MBSF10_10
Rigid bodies/nodes
9
100
53381
Model DOFs
54
600
160143
Computing time (s)
0.00051
0.042
18.41Figure11 presents the MAC matrices that show mode shape similarity of MBSF3_3 vs. FEA modes and MBSF10_10 vs. FEA modes, respectively. In contrast with that in Figure 11(a), the MAC matrix in Figure 11(b) has high diagonal values that are close to 1, excepting for the 4th and 5th modes, which indicates that the mode shapes correlated well between the two mode sets from the MBSF10_10 and the FE model. It is worth noting that mode 4 and 5 computed by the FE model actually have a repeated eigenvalue for the structural symmetry; thus, the related eigenvectors are not unique and hard to be matched.Figure 11
MAC matrix showing mode shape similarity between modes from MBSF models and the FE model: (a) MBSF3_3 vs. FEA and (b) MBSF10_10 vs. FEA.
(a)
(b)Figure11 also implies the correspondence between the MBSFs and FEA modes, which connects the similar mode-pairs together. Based on this correlation, frequency relative errors of the first 10 nonzero modes with respect to the FEA results are calculated and listed in Table 3. As can be seen that the maximum absolute value of the relative errors between the MBSF10_10 and the FEA modes is less than 5%, the average value is 3.4%.Table 3
Relative errors of mode frequencies with respect to FEA results.
Mode
Frequency (Hz)
Absolute value of relative error (%)
FEA
MBSF3_3
MBSF10_10
MBSF3_3
MBSF10_10
1
377.4
388.7
394.7
3.0
4.6
2
556.1
568.6
588.0
2.2
5.7
3
674.2
618.0
644.9
8.3
4.4
4
964.8
823.1
976.5
14.7
1.2
5
964.9
943.3
990.5
2.2
2.7
6
1662.6
1250.3
1628.0
24.8
2.1
7
1662.6
1310.5
1629.6
21.2
2.0
8
1735.5
1303.0
1717.4
24.9
1.0
9
1911.3
1735.8
2001.0
9.2
4.7
10
2049.9
1520.4
1933.7
25.8
5.7To compare the vibration behaviors of the two types of models comprehensively, the force-acceleration transfer curves fromz-direction of point P1 to z-direction of point P25 are plotted in Figure 12. It can be observed that the transfer function curve calculated by the MBSF10_10 (marked with solid magenta line) is in good agreement with the curve calculated by the FE model (marked with solid blue line). By comparison, a conclusion can be drawn that the vibration behavior of MBSF constructed by the RMMS, while divided reasonably, is consistent with the FE model.Figure 12
Comparison of force-acceleration transfer functions between 200 and 3000 Hz: P1-z to P25-z.


## 5.1.4. Comparison of MBSMs Results and Experiment Results

This section further compares the experimental modal results with the simulated ones computed by MBSMs, whose parameters are identified from experimentally modified FEA results. The model information of the MBSMs is listed in Table4.Table 4
Model scales of the MBSMs.
Items
MBSM3_3
MBSM6_6
MBSM8_8
MBSM10_10
Rigid bodies
9
36
64
100
Model DOFs
54
216
384
600Figure13 presents the stabilization diagram of the experimental normal modes in the LMS test. The mode frequencies and related shapes can be obtained in a laboratory. The diagram shows that the measured (MEA) results have good consistency and high reliability for comparison.Figure 13
Stabilization diagram of the experimental results.Figure14 presents the MAC matrices showing the mode shape correlation between the modes of the experiments and that of MBSMs with different division schemes, ranging from 3×3 to 10×10. It can be seen that mode shapes of the MBSM with an insufficient dividing scheme, such as MBSM3_3, are poorly correlated to the experimental ones because its MAC matrix has many high off-diagonal values, whereas good consistencies are achieved under reasonable division since MBSM6_6.Figure 14
MAC matrix showing mode shape similarity between modes from MBSM models and experimental results: (a) MBSM3_3 vs. MEA, (b) MBSM6_6 vs. MEA, (c) MBSM8_8 vs. MEA, and (d) MBSM10_10 vs. MEA.
(a)
(b)
(c)
(d)According to the correspondence between the MBSMs and FEA modes, as implied in MAC graphs of Figure14, their frequency relative errors in the first 10 nonzero modes are calculated and plotted in Figure 15. It can be seen that the maximum relative errors of the MBSMs with divisions from 6×6 to 10×10, with respect to the MEA results, are all less than 10%, and most of them are small than 5%.Figure 15
Mode frequency relative errors with respect to experimental results.In order to intuitively compare the vibration behaviors obtained by the simulation models and the experimental data, the force-acceleration transfer curves from different input channels to different out channels are plotted in Figures16 and 17. The results show that the transfer curves computed by the MBSMs are in good agreement with the experimental ones. It can be concluded that the vibration behaviors of MBSFs, which are established by the RMMS with experimental modification in parameter identification, are consistent with the measured results. Besides, the higher the division density, the better the consistency. However, the increase of rigid elements obviously amplifies the computational burden.Figure 16
Comparison of force-acceleration transfer functions between 200 and 1500 Hz: P1-z to P25-z.Figure 17
Comparison of force-acceleration transfer functions between 200 and 1500 Hz: P9-z to P21-z.


## 5.2. Case II: Application in a Lithography System

An additional application in concern with the dynamic modeling of lithography, which is key to the system design, was presented. The employed lithography, as shown in Figure18(a), mainly consists of an illuminator system, a reticle system, a wafer system, a projection system, precision measurement systems, vibration isolation systems, and frames. Both the reticle system and the wafer system are composed of a coarse stage and a fine stage. Among them, the coarse stage is to realize the micron-level high-acceleration motion, while the fine stage is to achieve the nanolevel ultraprecision motion. The illuminator system, reticle system, wafer system, projection system, and precision measurement system are all installed on the main frame. Also, three vibration isolators are used to connect the main frame with the base frame (i.e., the external world), so as to guarantee an ultrastable environment. In addition, there exist many flexible connections to further isolate the vibration, such as between the main frame and the projection lens, between the main frame and the supporting structures of the two motion stages.Figure 18
Structure of the lithography system: (a) schematic diagram and (b) rigid multibody model.
(a)
(b)In modeling of key structures/subsystems such as the illuminator system, the reticle, and wafer stages, the proposed parameter-identification method and RMMS were applied to establish the equivalent MBSs with parameters identified from modal results of related FE models, as shown in Figure19. Finally, a rigid multibody model of the lithography, which has 75 bodies and 127 DOFs, was constructed and presented in Figure 18(b). With this established model, the modal behavior was analyzed and compared with the experimental results.Figure 19
Dynamic models of key components/subsystems: (a) FEA models and (b) equivalent MBSs.The simulation and experimental modal results within 200 Hz are compared in Table5; it can be observed that the maximum relative error of mode frequency is about 6% while the mode shapes are matched. In addition, Figure 20 compares the simulation and experimental mode frequencies within 500 Hz. It can be summarized that vibration behavior computed by the final established rigid multibody model shows good agreement with the results from measured data, and only small errors exist between measured and the corresponding computational mode orders, especially within 0∼200 Hz, which is a mainly concerned frequency band in the design of lithography. It should also be noted that some modes are not caught in experiments due to the arrangement of sensors and limits of excitations.Table 5
Comparison of simulation and experimental modal results within 200 Hz.
Mode
Mode freq. (Hz)
Error (%)
Mode shape
Simulation
Experiment
Object
Mode description
1
1.49
1.5
0.4
Main frame
Translation alongY
2
1.5
1.51
0.5
Main frame
Translation alongX
3
1.99
2
0.1
Main frame
Translation alongZ
4
2.59
2.61
0.8
Main frame
Rotation aroundX
5
2.67
2.62
1.9
Main frame
Rotation aroundY
6
2.98
3
0.1
Main frame
Rotation aroundZ
7
37.06
39.01
5
Support plate of reticle stage
Torsional vibration
8
39.39
39.03
1
Support plate of reticle stage
Torsional vibration
9
39.47
39.0
1.3
Support plate of reticle stage
Torsional vibration
10
49.7
50
0.6
Wafer stage
Translation alongX
11
69.19
70
1.2
Wafer stage
Translation alongX + rotation around Z
12
76.87
78
1.5
Wafer stage
Rotation aroundX
13
80.33
80
0.5
Reticle stage
Translation alongY
14
82.52
81
1.9
Wafer stage
Rotation aroundZ + translation along X
15
101.9
98
4
Reticle stage
Translation alongX
16
112.9
120
6
Wafer stage
Rotation aroundZ
17
124
128
3.2
Reticle stage
Rotation aroundX
18
127.3
134
5
Lens + reticle stage
—
19
128.4
—
—
Reticle stage
Translation alongZ
20
147.8
150
1.4
Reticle stage
Rotation aroundX
21
171.2
173
1
Lithography
—
22
176.8
182
2.7
Reticle stage
Rotation aroundZ
23
190.8
195
2.2
Reticle stage
Rotation aroundYFigure 20
Comparison of simulation and experimental modal frequencies within 500 Hz.


## 6. Discussion

The abovementioned experiments and application strongly illustrate the feasibility and correctness of the proposed parameter identification method and the RMMS to predict the dynamics of flexible structures by establishing equivalent rigid MBSs with less DOFs and explicit physical meaning.The simulation results in Case I show the RMMS can be applied in static analysis since the deformations fit well with the FEA results. Meanwhile, the further comparisons show that dynamics of the established MBSs not only fit well to the FE model without modification in parameter identification but also fit well to the experimental ones while the substituted parameters are identified from experimentally modified FEA results. In addition, the results fit better as the division density grows within a certain range. Nevertheless, a trade-off is necessary between the high accuracy and low model DOFs when determining the division scheme. Generally, the accuracy of FE models can be improved by experimental modification as well, but still yielding high order models. Therefore, it can be concluded that the RMMS has the advantage of low DOFs in comparison with the FE method, whereas the explicit physical meaning derives benefits of good reusability and extensibility, in comparison with the matrix-based MOR methods.Application in a lithography system also shows that the parameter identification method and the RMMS is suitable for a large-scale mechanical system since only the component- or subsystem-level FE models of few less-rigid structures are required rather than the whole system-level FE model. While for the case of common structures that can be modeled by FE method with high reliability, the modification may be omitted, for the case of complex structures, a more recommended way is to update the target eigenpairs by the measured data if possible. Once the final equivalent rigid MBS is established, it can overcome the low computational efficiency and poor real-time problems to serve for dynamic analysis, iterative optimization, and even the control in ultraprecision systems.


## 7. Conclusions

In this paper, an RMMS, which directly discretizes complex flexible structures into RFEs connected by SDEs, was proposed. The core parameter identification problem in the RMMS was expounded in detail and resolved from the perspective of the inverse vibration problem. Through the proof of the equivalence of SDEs at arbitrary position and orientation between two rigid bodies, the nonlinear inverse problem inRn×n matrix space was decoupled and linearized into solving the incompatible linear equations in Rn2 vector space based on the vectorization operator and Kronecker product. The optimal dynamic parameters were eventually identified by calculating the Moore-Penrose generalized inverse. Finally, the agreements of the experiment results on a planar motor mover with simulated ones by the RMMS strongly validated the feasibility and correctness of the proposed parameter identification method and the RMMS in accurately predicting the dynamic behaviors, while with few DOFs and explicit physical meaning. The application in a lithography system indicates that the RMMS is applicable for modeling the flexibility of weak-rigid structures in the dynamic analysis of large-scale mechanical systems, especially in dynamic design and iterative optimization.



---

*Source: 8329395-2020-01-28.xml*