# Modeling and Analysis of a Multi-Degree-of-Freedom Micro-Vibration Simulator

Authors: Xiaoming Wang; Zhenbang Xu; Shuai He; Qingwen Wu; Hang Li; Yaping Zhao; Lei He
Journal: Shock and Vibration
 (2017)

---

## Abstract

To reproduce the disturbance forces and moments generated by the reaction/momentum wheel assembly, a multi-degree-of-freedom micro-vibration simulator is proposed. This can be used in the ground vibration experiments of an optical payload replacing the real action/momentum wheel assembly. First, the detailed structure of the micro-vibration simulator is introduced. Then, the complete system kinematic and dynamic models of the micro-vibration simulator are derived. In addition, the disturbance forces and moments produced by the micro-vibration simulator are calculated. Finally, the normal mode analysis and a cosimulation are adopted to verify the validity of this method. The analysis and simulation results show that the micro-vibration simulator can exactly reproduce the disturbance forces and moments with different amplitudes and different frequency ranges.


---

## Body

## 1. Introduction

With the development of space technology, space optical remote sensors with large aperture and high-resolution imaging are playing an increasing role in many fields of military, science, astronomy, and civil use. However, micro-vibrations produced by the operation of on-board equipment, such as reaction/momentum wheel assemblies (R/MWAs), thrusters, electric motors, and cryogenic coolers, can cause significant effects in the image quality degradation [1].Among the various micro-vibration disturbance sources, that produced by R/MWAs is generally regarded as the largest [2–4]. There have been numerous studies on R/MWA disturbance modeling, as well as the characteristics of micro-vibration caused by R/MWAs [2, 3, 5]. The general approach to reducing the disturbance effect on the image quality of an optical payload is to insert one or several isolators in the transmission path between the disturbance source and the payload [6]. To develop the appropriate isolators and ensure mission success, a significant amount of analysis and ground experiments needs to be devoted to validating the isolation effect of the isolators and the image performance of optical payload before launch. During the ground experiments, the effects of R/MWA induced disturbance on the optical payload are of prime concern. Zhou et al. [7, 8] and Kamesh et al. [9, 10] all used a real momentum wheel assembly (MWA) to measure the micro-vibrations produced by the MWA with and without the proposed soft suspension system. However, it is not usual to employ all the flight R/MWAs to conduct the ground validating experiments because of scheduling issues or product assurance activities. Therefore, developing a micro-vibration simulator, which can replace the real flight R/MWAs, is considered to be an important adjunct to the development processes for space missions.Hostens et al. [11, 12] have proposed a six-degree-of-freedom (DOF) vibration simulator driven by hydraulic, which can be used to test the vibration of mobile machinery. The vibration simulator adopted Gough-Stewart platform, which greatly improved the structural stability. Parallel manipulator such as the Gough-Stewart platform (GSP) has been widely used in many different application fields for its remarkable performance (e.g., high maneuverability, precision, and high load/weight capacity). Dynamic analysis of the parallel manipulator plays a vital role in the design and control of such manipulators. Over the past decades, numerous research results have been reported on the kinematics of the GSP and relatively fewer results on the dynamic of GSP [13]. Oftadeh et al. [14] presented the explicit dynamics formulation for the GSP and utilized the Lagrange method to verify the resulting dynamics equations. Dasgupta and Mruthyunjaya [15] derived an inverse dynamic formulation by the Newton-Euler approach for the GSP, with the frictional forces occurring in the joints; meanwhile the mass of inertia of the pods was taken into consideration in their study. Pedrammehr et al. [16] produced the improved dynamic equations of the GSP, which took into account the rotational degree of the pods around the axial direction. Similar efforts can be found in [17–20]. Generally in the literature, many researches focused on the influence of the inertia forces, Coriolis forces, and the frictional forces on the dynamics control of the GSP. However, researches on the disturbance forces and moments acting on the mounting surface of the GSP are still rare. Moreover, the calculation of disturbance forces and moments is critical for the GSP based micro-vibration simulator.There are two kinds of simulation methods of micro-vibration. One is outputting micro-vibration through the moving platform. The first method has been studied by many scholars. For example, Hostens et al. [11, 12] designed a vibration simulator which can simulate 6-dimensional vibrations. The other is outputting disturbance forces and moments through the base platform. A few scholars have studied this method. Park et al. [21] developed a micro-vibration emulator to test the jitter in spacecraft. However, this emulator can only produce three disturbance forces which are mutually perpendicular. Therefore, developing a micro-vibration simulator which can reproduce 6-dimensional disturbance forces and moments is novel.In this study, a multi-degree-of-freedom micro-vibration simulator (MMVS) based on the GSP, which can exactly reproduce the 6-dimensional disturbance forces and moments, is presented. The dynamic models of the MMVS, which consider the effects of actuator inertia and eccentric load, are derived using the Newton-Euler method and Lagrange approach. The formulation incorporates all the elastic, inertia, Coriolis, centrifugal, and external forces; furthermore, the disturbance forces and moments produced by the MMVS are also analyzed. Based on the derived dynamic models, a disturbance planning method is proposed to obtain the target disturbance forces and moments. The cosimulation using ADAMS and MATLAB®/Simulink is adopted to verify the validity of the dynamic models and the feasibility of MMVS.This paper is organized as follows: in Section2, we describe the structure of the MMVS. Section 3 presents the dynamics model including details of the calculation of actuator forces and planning method of the disturbance forces and moments. Section 4 presents the normal mode analysis and cosimulation verification. Concluding remarks on the micro-vibration simulator are presented in Section 5.


## 2. Micro-Vibration Simulator Description

The detailed structure of the MMVS is depicted in Figure1. It consists of a moving platform, base plane, and six actuators. Each actuator, driven by a voice coil motor, is connected to the base plane by a universal joint at one end and the moving platform by a gimbal, which has three degrees-of-freedom, at the other end. Figure 2 shows a cross-sectional view of the actuator with a permanent magnet attached to the envelope. The two membranes perform the function of spring only for axial compliance. A connecting rod attached to the central point of the upper membrane connects the lower membrane and the voice coil. The one active degree-of-freedom (DOF) on the actuator is the axial movement of the coil.Figure 1
Virtual prototype of the vibration simulator.Figure 2
Cross-sectional view of actuator.The working principle of this proposed MMVS is as follows: the control signal produced by the upper computer is amplified by the power amplifier to drive the voice coil motors. The reciprocating motion of the actuators causes the moving platform to move together. As a result, the inertial forces are generated. According to Newton’s Third Law, the inertial forces produced by the moving parts will be transmitted to the mounting surface through the universal joints, which leads to the generation of disturbance forces and moments. The object of this study is to calculate the actuator forces, which are used to generate the target disturbance forces and moments.


## 3. Dynamics Model

The object of this study is to develop a MMVS which can reproduce the desired disturbance forces and moments. So we need to calculate the actuator forces which can be used to derive the MMVS to obtain the desired disturbance forces and moments. First, the dynamic relations between the actuator forces and the generalized acceleration, velocities, and displacements are derived. It is noted that the assumption is adopted that the generalized velocities and generalized displacements of the moving platform are the time integral and time double integral of the generalized acceleration, respectively. Then we calculate the disturbance forces and moments of the MMVS acting on the mounting surface when the generalized acceleration of the moving platform is given. Finally the dynamic relations between the actuator forces and the desired disturbance forces and moments are carried out. The planning method of disturbance forces and moments is depicted in Figure3.Figure 3
Block diagram of the planning method of disturbance forces and moments.


### 3.1. Kinematics

The scheme of the vibration simulator is shown in Figure4, in which the coordinate system {P} is the body frame fixed to the geometric center of the joints of the moving platform, while the {B} coordinate system is the base frame attached to the geometric center of the joints of the base plane. The {g} coordinate system is the inertial frame fixed to the geometric center of the bottom of the base plane, and its orientation is identical with frame {B}. The linear motions are denoted as surge (x), sway (y), and heave (z) along theXB–YB –ZB axes of the base frame, and the angular motions roll (γ), pitch (β), and yaw (α) areX-Y-Z fixed angles. The upper joint points and the lower joint points are denoted with pPi in frame {P} and pBi in frame {B}, respectively. RP and RB describe the radii of the payload and base platforms. The angle between P6 and P1 is denoted by φ. The angle between B6 and B1 is denoted by θ, which is illustrated in Figure 4(b).Figure 4
A schematic view of the vibration simulator: (a) isometric view and (b) vertical view.
(a)
(b)


#### 3.1.1. Velocity Analysis

In the base frame, the kinematic equations of theith actuator can be described by(1)lBi=tB+pBi-bBi=tB+RPB·pPi-bBi,where the variable i defines the ith actuator variables, indicating that in general the equations are applicable to any actuator. lBi is the length vector of the actuator with respect to the base frame {B}, tB is the position vector of the body frame, tB=xyzT, and RPB is the rotation matrix of the transformation from the body frame {P} to base frame {B}. The rotation matrix is given by(2)RPB=cos⁡αcos⁡βcos⁡αsin⁡βsin⁡γ-sin⁡αcos⁡γcos⁡αsin⁡βcos⁡γ+sin⁡αsin⁡γsin⁡αcos⁡βsin⁡αsin⁡βsin⁡γ+cos⁡αcos⁡γcos⁡αsin⁡βcos⁡γ-cos⁡αsin⁡γ-sin⁡βcos⁡βsin⁡γcos⁡βcos⁡γ.Taking the derivative of (1) with respect to time, the velocity mapping function can be obtained and is given by(3)q˙pi=l˙Bi=t˙B+ω×RPBpPi=E3RPBp~PiTRPBTq˙,where q˙pi is the velocity of the joint Pi, t˙B and ω are the translational velocity and angular velocity of the moving platform, respectively, t˙B=x˙y˙z˙T, and ω=γ˙β˙α˙T, E3 is a unit 3 × 3 matrix, p~Pi is the skew symmetry matrix of pPi, q˙ is the general velocity of the moving platform, and q˙=t˙BTωTT.LettingJpiq=E3RPBp~PiTRPBT, (3) can be rewritten as(4)q˙pi=Jpiq·q˙,where Jpiq denotes the Jacobian matrix relating the general velocity to the velocity of the upper joint pi.The length of actuator can be described as(5)li=lBiT·lBi.Taking the time derivative of (5), the sliding velocity of the actuator can be obtained:(6)l˙i=lBniT·t˙B+lBniT·ω×RPB·pPi,where lBni is the unit vector of lBi.Equation (6) can be rewritten as(7)l˙i=lBniT,RPB·pPi×lBniT·t˙Bω.Thus (8) may be obtained as(8)L˙=J·q˙,where L˙ denotes actuator sliding velocities, L˙=l˙1,l˙2,l˙3,l˙4,l˙5,l˙6T, and J is the actuator Jacobian matrix expressing mappings from the general velocity to the actuator sliding velocities, which is given by(9)J=lBn1T,RPB·pP1×lBn1TlBn2T,RPB·pP2×lBn2TlBn3T,RPB·pP3×lBn3TlBn4T,RPB·pP4×lBn4TlBn5T,RPB·pP5×lBn5TlBn6T,RPB·pP6×lBn6T.According to their physical meaning, the velocities of the actuator can also be described in terms of the sliding and angular velocities of the actuator underframe{B} as(10)l˙Bi=l˙i·lBni+ωai×lBi,where ωai is the angular velocity of the actuator.Taking the cross product of the above equation withlBni and considering the assumption that no rotation is allowed about the actuator axis (i.e., ωaiT·lBi=0) yield(11)ωai=lBni×l˙Bili=lBni×q˙pili.


#### 3.1.2. Acceleration Analysis

For convenience of analysis, the actuator is derived into two parts: the upper leg and the lower leg as shown in Figure5.Figure 5
Diagram of the actuator.Taking the derivative of both sides of (11) with respect to time and taking into account the fact that ωaiT·lBni=0 and εaiT·lBni=0, the angular acceleration of the actuator can be obtained as(12)εai=lBni×q¨pi-2ωai·l˙ili.The translational acceleration of the upper jointPi can be obtained by taking the derivative of (3) with respect to time as described by(13)q¨pi=t¨B+ω˙×RPB·pPi+ω×ω×RPB·pPi.According to their physical meaning, the velocities of the upper jointPi can be described in terms of the velocity of the center of mass of the upper leg of the underframe {B} as(14)q˙pi=vBrci+ωai×rrci·lBni.Rearranging (14) yields(15)vBrci=q˙pi-ωai×rrci·lBni.Substituting (11) into (15), the velocity of the upper leg centroid is given by(16)vBrci=q˙pi+lBni×lBni×q˙pirrcili.LetJrcipi=E3+rrci/li·l~Bni2, and (16) can be rewritten as(17)vBrci=Jrcipiq˙pi,where Jrcipi denotes the Jacobian matrix relating the velocity of the upper joint pi to the velocity of the upper leg centroid.In the following section, to obtain the acceleration of the upper leg centroid, (18) can be used, which is given by(18)dl~BniT·lBnidt=-l~BniTl~Bniq˙pilBniT+lBniq˙piTl~BniTl~Bnili.Taking the derivative with respect to time, of (17), and considering (18), the acceleration of the upper leg centroid is given by(19)aBrci=Jrcipi·q¨pi+2l~BniTl~Bniq˙pilBniTq˙pi+lBniq˙piTl~BniTl~Bniq˙pirrcili2.The same procedure may be easily adopted to obtain the velocity and acceleration of the lower leg, which are given by(20)vBtci=rtci·ωai×lBni=Jtcipi·q˙pi,(21)aBtci=Jtcipi·q¨pi-2l~BniTl~Bniq˙pilBniTq˙pi+lBniq˙piTl~BniTl~Bniq˙pirtcili2,where Jtcipi denotes the Jacobian matrix relating the velocity of the upper joint pi to the velocity of the lower leg centroid, Jtcipi=l~BniTl~Bnirtci/li.When the center of mass of the moving platform does not coincide with the origin of the frame{P}, the centrifugal and Coriolis forces will manifest themselves. Therefore the centrifugal and Coriolis terms cannot be ignored when we derive explicit equations for the dynamics of the moving platform. Let pPc denote the position vector of the centroid of the moving platform under body frame {P}, which can be shown by(22)pBc=tB+RPB·pPc.Taking the derivative of both sides of (22) with respect to time, the acceleration of the centroid of the moving platform is given by(23)p¨Bc=t¨B+RPB·p~PcT·RPBT·ω˙+ω~2·RPB·pPc.So far, all the kinematic parameters required for the inverse dynamics are derived. In the next section, attention will be paid to the deduction of dynamic equations of the MMVS.


### 3.2. Actuator Forces

In this subsection, the Lagrange approach and Newton-Euler method are adopted to develop the dynamic equations of the MMVS. Because of the function of the membrane of the actuators, the effect of gravity on the moving platform and the actuators can be compensated by the elastic forces of the membrane; that is, the moving platform will move to a new equilibrium position; therefore the effect of gravity on the micro-vibration simulator may be ignored.


#### 3.2.1. Analysis of a Single Actuator

For the dynamic modeling of the actuator, only the actuator forcesfi and the interaction forces fpi between the upper joints are covered in this subsection. Figure 6 shows a schematic of the actuator. As shown in Figure 6, frame s^i is the body frame of the upper leg which is attached to the center of mass of upper leg with s^1i along theith actuator. Here the over hat indicates a unit of length. The d^i coordinate system is located to the ith actuator at the mass center of the lower leg with d^1i along theith leg.Figure 6
Schematic of the actuator.The Lagrange equation of the actuator can be written as in(24)d∂Ti/∂q˙pidt-∂Ti∂qpi=Qi,where Qi is the generalized force component associated with the generalized coordinate qpi and T is the kinematic energy of the ith actuator, given by(25)Ti=vBrciTmrciBvrci2+vBtciTmtciBvtci2+ωaiTIBrci+IBtciωai2,where IBrci and IBtci are the moments of inertia of the upper leg and lower leg with respect to frame {B}, respectively.The moments of inertiaIBrci and IBtci are discussed in what follows. Assuming that the actuator is axisymmetric and the principal axes of inertia of the upper leg are consistent with the corresponding axes of frame s^i, the principal moments of inertia are denoted by Is^1s^1, Is^2s^2, and Is^3s^3, and the following expression can be obtained:(26)IBrci=RBs^iIs^irciRBs^iT=RBs^i·diag⁡Is^1s^1iIs^2s^2iIs^3s^3i·RBs^iT=Is^1s^1is^1is^1iT+Is^2s^2is^2is^2iT+Is^3s^3is^3is^3iT,where RBs^i is the rotation matrix of the transformation from frame s^i to base frame {B}, given by(27)RBs^i=s^1iT·i^s^2iT·i^s^3iT·i^s^1iT·j^s^2iT·j^s^3iT·j^s^1iT·k^s^2iT·k^s^3iT·k^,where i^=100T, j^=010T, and k^=001T are three unit vectors, which are used to present the base frame {B}.Because the actuator is axisymmetric, the principal moment of inertiaIs^2s^2i is equal to Is^3s^3i. Letting Iri denote Is^2s^2i, (26) can be rewritten by(28)IBrci=Is^1s^1is^1is^1iT+IriE3-s^1is^1iT.Similarly, the moment of inertiaIBtci can be described by(29)IBtci=Id^1d^1id^1id^1iT+ItiE3-d^1id^1iT,where Id^1d^1i is the principal moment of inertia of the lower leg with respect to the axis d^1i of frame d^i and Iti presents the principal moment of inertia of the lower leg with respect to the axis d^2i (or the axis d^3i) of frame d^i.Substituting (17), (20), (28), and (29) into (25) yields(30)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12ωaiTIs^1s^1is^1is^1iT+IriE3-s^1is^1iT+Id^1d^1id^1id^1iT+ItiE3-d^1id^1iTωai.From Figure6, it can be known that s^1i, d^1i, and lni are equal to each other. Thus, (30) can be rewritten as(31)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12Is^1s^1i+Id^1d^1i-Iri-ItiωaiT·lBni2+Iri+ItiωaiT·ωai.Considering the assumption that no rotation is allowed about the actuator axis (i.e.,ωaiT·lBni=0) and substituting (11) into (31), the following expression can be obtained: (32)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12Iri+ItiωaiT·ωai=12q˙piTJrcipimrci+Jtcipimtci+l~BniT·l~BniIri+Itili2q˙pi=12q˙piTMlegiq˙pi,where Mlegi=Jrcipimrci+Jtcipimtci+l~BniT·l~BniIri+Iti/li2.The first and second terms in (24) can be written as(33)d∂T/∂q˙pidt=Mlegiq¨pi+q˙pi2lBniq˙piTl~BniTl~Bnimrcirrcili2+lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT-2lBniTq˙piTl~BniTl~Bni+lBniq˙piTl~BniTl~Bni+l~BniTl~Bniq˙pilBniTmrcirrci2+mtcirtci2+Iri+Itili3,∂T∂qpi=lBniq˙piTq˙pi+2lBniTq˙piq˙pi-3lBniq˙piTlBnilBniTq˙pimrcirrcili2-lBniq˙piTq˙pi+q˙piq˙piTlBni-2lBnilBniTq˙piq˙piTlBnimrcirrci2+mtcirtci2+Iri+Itili3.The right side of (24) can be described as(34)Qi=lBnifi+fpi.Substituting (25), (33), and (34) into (24) yields(35)fpi=Mlegiq¨pi+Caiq˙pi-lBnifi=MlegiJpiqq¨+CaiJpiqq˙+Mlegiω~2RPBpPi-lBnifi,where Cai is a 3×3 matrix, which is given by(36)Cai=mrcirrcili2·lBniq˙piTl~BniTl~Bni+lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT-mrcirrci2+mtcirtci2+Iri+Itili3·lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT.


#### 3.2.2. The Completed Dynamic Equations

Assuming that the axial stiffness coefficients of the actuators are equal and lettingk denote the axial stiffness coefficient, then for theith actuator (37)fei=-k·Δli,where fei is the force needed to cause a Δl change in theith actuator length. Assembling the equations for all the actuators and considering (8), (37) becomes(38)fe=-k·J·Δq,where fe is the general elastic force, fe=fe1fe2fe3fe4fe5fe6T.Similarly, the general damping forcefc can be described as (39)fc=-c·J·q˙,where fc is a 6×1 matrix and c is the damping coefficient of the actuator.Based on the virtual work principle, the dynamic equation of the moving platform under the external forces of the loads can be described as(40)ΓT·δq=feT+fcT+FT·δL,where Γ=τ1,τ2,τ3,τ4,τ5,τ6T is the applied generalized force exerted on the moving platform, δq=δx,δy,δz,δγ,δβ,δαT is virtual displacement of the moving platform, F=f1,f2,f3,f4,f5,f6T describes the actuator forces, and δL=δl1,δl2,δl3,δl4,δl5,δl6T gives the virtual displacement of the actuators.(41)δL=J·δq.Substituting (41) into (40), the external force Γ acting on the moving platform can be written as(42)Γ=JTfe+fc+F.Substituting (38) and (39) into (42) yields(43)Γ=-K·Δq-C·q˙+JTF,where K=kJTJ is the matrix for generalized stiffness and C=cJTJ is the matrix for generalized damping.The dynamic model for upper platform as a rigid body ignoring the inertia of the actuators can be derived using the Newton-Euler method [14]:(44)m·p¨Bc=F3×1,RPB·IP·RPBT·ω˙+ω~·RPB·IP·RPBT·ω+m·RPB·pPc×t¨=M3×1,where m is the payload mass, IP is the inertia matrix with respect to frame {P}, F3×1=τ1,τ2,τ3T and M3×1=τ4,τ5,τ6T.Substituting (43) into (44) yields(45)mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBTq¨+000ω~RPB·IP·RPBT+c·JTJq˙+mE30ω~2·RPB·p~Pc+k·JT·J·Δq=JTF.Considering the inertia of the actuators, (45) can be written as(46)mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBTq¨+000ω~RPB·IP·RPBT+c·JTJq˙+mE30ω~2·RPB·p~Pc+k·JT·J·Δq=-HpFP.Expressions forHp and Fp are given by(47)Hp=E3E3E3E3E3E3RPBp~P1RPBTRPBp~P2RPBTRPBp~P3RPBTRPBp~P4RPBTRPBp~P5RPBTRPBp~P6RPBT,FP=fp1Tfp2Tfp3Tfp4Tfp5Tfp6TT.Substituting (35) and (47) into (46) produces the complete dynamic equation of the vibration simulator, which can be described as(48)Mqq¨+Cq,q˙q˙+K·Δq=JTF,where M(q) is a 6×6 mass matrix, Cq,q˙ is a 6×6 matrix of the centrifugal and Coriolis force terms, K is a 6×6 matrix of the generalized stiffness, and F is a 6×1 vector representing actuator forces. They are given by(49)Mq=mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBT+∑i=16JpiqTMlegiJpiq,Cq,q˙q˙=000ω~RPB·IP·RPBT+c·JTJq˙+∑i=16JpiiTCaiJpiiq˙+mE30ω~2RPB·pPc+∑i=16JpiiTMlegiω~2RPB·pPi.Premultiplying both sides of (48) with J-T, the actuator forces are given by(50)F=J-TMqq¨+Cq,q˙q˙+K·Δq.


### 3.3. Calculation of Disturbance Forces and Moments

As we know, the object of this study is to develop a micro-vibration simulator which can reproduce desired disturbance forces and moments. So this subsection will calculate the disturbance forces and moments when the actuator forces are applied to the MMVS.The calculation of the disturbance forces and moments is based on the Newton-Euler method. For the convenience of analysis, a frame of reference is attached to the actuator with its origin at the lower joint points,x-axis (hi) along the actuator, y-axis (yi) along the rotation axis (axis fixed to the actuator) of the universal joint, andz-axis (ki) perpendicular to the x- and y-axes according to the right hand rule. Figure 1 shows the frame. Considering the rotational equilibrium of the entire actuator, Euler’s equation for the actuator is obtained as(51)IBrci+IBtciεai+ωai×IBrci+IBtciωai+mrcirrci×aBrci+mtcirtci×aBtci=lBi×fpi+Muihi,where rrci=(li-rrci)lBni, rtci=rtcilBni, and Mui is the magnitude of the constraint moment at the universal joint acting about the hi axis.Taking cross products of both sides of (51) with hi, the scalar unknown Mui is given by (52)Mui=lBniTCilBniThi,where Ci is a 3×1 matrix, which is given by(53)Ci=IBrci+IBtciεai+ωai×IBrci+IBtciωai+mrcirrci×aBrci+mtcirtci×aBtci.Ignoring the damping forces of the actuator and considering the upper leg of the actuator, Newton’s equation can be given by(54)fi+lBniTfpi-k·Δli-mrcilBniT·aBrci=0.Solving (54), the constraint force at the spherical joint can be rewritten as(55)fpi=XilBni+Di.ExpressionsXi and Di are given by(56)Xi=mrcilBniT·aBrci+k·Δli-fi,Di=Ci×lBni-Muihi×lBnili.Considering the translational equilibrium, Newton’s equation for the entire actuator can be written as(57)mrciBarci+mtciBatci-fpi=fui,where fui is the constraint force at the universal joint.Substituting (55) into (57), the constraint force fui is given by(58)fui=mrciBarci+mtciBatci-XilBni+Di.Consequently, the disturbance forces acting on the center of the bottom of the base plane, which are produced by the simulator, can be written as(59)Fm=∑i=16fui.Similarly, the disturbance moments acting on the center bottom of the base plane, produced by the simulator, which are shown in Figure6, can be written as(60)Mm=∑i=16rOBi×fui+Muihi,where rOBi is the vector from the origin Og to the lower joint point Bi.


### 3.4. Planning of Disturbance Forces and Moments

In the last subsection, the relationships between the disturbance forces and moments and the actuator forces are derived. This subsection will deduce the relationship between the deduced forces and moments and the generalized acceleration of the moving platform and then plan the disturbance forces and moments.Because micro-vibrations have the characteristics of a wide frequency range and small amplitudes, a few approximations can be adopted to plan the disturbance forces and moments. For example, the third term on the right side of (13) is a square term of angular velocity, so it can be ignored when the angular velocity of the moving platform is small. The simplified acceleration of the upper jointPi can be described as(61)p¨Bi≅t¨B+ω˙×RPB·pPi≅E3p~BiTt¨Bω˙=Tpiq¨,where Tpi=E3p~BiT.Similarly, ignoring the second terms of (19) and (21), respectively, the acceleration of the upper and lower legs is given by (62)aBrci=Jrcipi·q¨pi,aBtci=Jtcipi·q¨pi.Substituting (61) into (62) yields(63)aBrci=JrcipiE3p~BiTt¨Bω˙=JrcipiTpiq¨,aBtci=JtcipiE3p~BiTt¨Bω˙=JtcipiTpiq¨.The angular velocity of the actuator is very small, so it can also be ignored. The angular acceleration of the actuator can be described as(64)εai≅l~BniaBpili.After obtaining the simplified expressions presented as (61), (63), and (64), the following expressions are obtained:(65)Xi≅mrcilBniT·JrcipiTpiq¨+k·Δli-fi,Ci≅IBrci+IBtci·l~Bnili+mrcir~rciJrcipi+mtcir~tciJtcipiTpiq¨,Di≅-l~BniliIBrci+IBtci/li·l~Bni+mrcir~rciJrcipi+mtcir~tciJtcipiTpiq¨=Λiq¨,where Λi=-l~Bni/liIBrci+IBtci/li·l~Bni+mrcir~rciJrcipi+mtcir~tciJtcipiTpi.Therefore, the constraint forces at the spherical and universal joints can be rewritten as(66)fpi=XilBni+Λiq¨,fui=mrciJrcipi+mtciJtcipi-mrcilBnilBniTJrcipiTpi-Λi+lBnieiJ-TMqq¨=Θuiq¨,where ei=[0⋯1⋯0] is a 1×6 matrix, theith element is 1, and the other elements are 0, and Θui is 3×1 matrix, Θui=mrciJrcipi+mtciJtcipi-mrcilBnilBniTJrcipiTpi-Λi+lBnieiJ-TMq.The disturbance forces and moments acting on the center of the base plane’s bottom can be rewritten as(67)Fm=∑i=16fui=∑i=16Θuiq¨=ΘTq¨,Mm=∑i=16rOBi×fui+Muihi=∑i=16r~OBiΘuiq¨=ΘRq¨.Equation (67) can be written in a matrix form as(68)FmMm=ΘTΘRq¨.Taking both sides of (68) with ΘTTΘRT-T, the generalized acceleration of the moving platform can be given by(69)q¨=ΘTΘR-1FmMm.LetFmd and Mmd denote the desired disturbance forces and moments, respectively; the desired generalized acceleration of the moving platform can be described by(70)q¨d=ΘTΘR-1FmdMmd.Assuming that the generalized velocities and generalized displacements of the moving platform are, respectively, the time integral and time double integral of the generalized acceleration, the following equations can be obtained:(71)q˙d=∫q¨ddt,qd=∬q¨ddt.Substituting (70) and (71) into (50) yields the desired actuator forces, which are given by(72)Fd=J-TMqdq¨d+Cqd,q˙dq˙d+Kqdqd.So far, the relations between the actuator forces and the desired disturbance forces and moments are obtained.


## 3.1. Kinematics

The scheme of the vibration simulator is shown in Figure4, in which the coordinate system {P} is the body frame fixed to the geometric center of the joints of the moving platform, while the {B} coordinate system is the base frame attached to the geometric center of the joints of the base plane. The {g} coordinate system is the inertial frame fixed to the geometric center of the bottom of the base plane, and its orientation is identical with frame {B}. The linear motions are denoted as surge (x), sway (y), and heave (z) along theXB–YB –ZB axes of the base frame, and the angular motions roll (γ), pitch (β), and yaw (α) areX-Y-Z fixed angles. The upper joint points and the lower joint points are denoted with pPi in frame {P} and pBi in frame {B}, respectively. RP and RB describe the radii of the payload and base platforms. The angle between P6 and P1 is denoted by φ. The angle between B6 and B1 is denoted by θ, which is illustrated in Figure 4(b).Figure 4
A schematic view of the vibration simulator: (a) isometric view and (b) vertical view.
(a)
(b)


### 3.1.1. Velocity Analysis

In the base frame, the kinematic equations of theith actuator can be described by(1)lBi=tB+pBi-bBi=tB+RPB·pPi-bBi,where the variable i defines the ith actuator variables, indicating that in general the equations are applicable to any actuator. lBi is the length vector of the actuator with respect to the base frame {B}, tB is the position vector of the body frame, tB=xyzT, and RPB is the rotation matrix of the transformation from the body frame {P} to base frame {B}. The rotation matrix is given by(2)RPB=cos⁡αcos⁡βcos⁡αsin⁡βsin⁡γ-sin⁡αcos⁡γcos⁡αsin⁡βcos⁡γ+sin⁡αsin⁡γsin⁡αcos⁡βsin⁡αsin⁡βsin⁡γ+cos⁡αcos⁡γcos⁡αsin⁡βcos⁡γ-cos⁡αsin⁡γ-sin⁡βcos⁡βsin⁡γcos⁡βcos⁡γ.Taking the derivative of (1) with respect to time, the velocity mapping function can be obtained and is given by(3)q˙pi=l˙Bi=t˙B+ω×RPBpPi=E3RPBp~PiTRPBTq˙,where q˙pi is the velocity of the joint Pi, t˙B and ω are the translational velocity and angular velocity of the moving platform, respectively, t˙B=x˙y˙z˙T, and ω=γ˙β˙α˙T, E3 is a unit 3 × 3 matrix, p~Pi is the skew symmetry matrix of pPi, q˙ is the general velocity of the moving platform, and q˙=t˙BTωTT.LettingJpiq=E3RPBp~PiTRPBT, (3) can be rewritten as(4)q˙pi=Jpiq·q˙,where Jpiq denotes the Jacobian matrix relating the general velocity to the velocity of the upper joint pi.The length of actuator can be described as(5)li=lBiT·lBi.Taking the time derivative of (5), the sliding velocity of the actuator can be obtained:(6)l˙i=lBniT·t˙B+lBniT·ω×RPB·pPi,where lBni is the unit vector of lBi.Equation (6) can be rewritten as(7)l˙i=lBniT,RPB·pPi×lBniT·t˙Bω.Thus (8) may be obtained as(8)L˙=J·q˙,where L˙ denotes actuator sliding velocities, L˙=l˙1,l˙2,l˙3,l˙4,l˙5,l˙6T, and J is the actuator Jacobian matrix expressing mappings from the general velocity to the actuator sliding velocities, which is given by(9)J=lBn1T,RPB·pP1×lBn1TlBn2T,RPB·pP2×lBn2TlBn3T,RPB·pP3×lBn3TlBn4T,RPB·pP4×lBn4TlBn5T,RPB·pP5×lBn5TlBn6T,RPB·pP6×lBn6T.According to their physical meaning, the velocities of the actuator can also be described in terms of the sliding and angular velocities of the actuator underframe{B} as(10)l˙Bi=l˙i·lBni+ωai×lBi,where ωai is the angular velocity of the actuator.Taking the cross product of the above equation withlBni and considering the assumption that no rotation is allowed about the actuator axis (i.e., ωaiT·lBi=0) yield(11)ωai=lBni×l˙Bili=lBni×q˙pili.


### 3.1.2. Acceleration Analysis

For convenience of analysis, the actuator is derived into two parts: the upper leg and the lower leg as shown in Figure5.Figure 5
Diagram of the actuator.Taking the derivative of both sides of (11) with respect to time and taking into account the fact that ωaiT·lBni=0 and εaiT·lBni=0, the angular acceleration of the actuator can be obtained as(12)εai=lBni×q¨pi-2ωai·l˙ili.The translational acceleration of the upper jointPi can be obtained by taking the derivative of (3) with respect to time as described by(13)q¨pi=t¨B+ω˙×RPB·pPi+ω×ω×RPB·pPi.According to their physical meaning, the velocities of the upper jointPi can be described in terms of the velocity of the center of mass of the upper leg of the underframe {B} as(14)q˙pi=vBrci+ωai×rrci·lBni.Rearranging (14) yields(15)vBrci=q˙pi-ωai×rrci·lBni.Substituting (11) into (15), the velocity of the upper leg centroid is given by(16)vBrci=q˙pi+lBni×lBni×q˙pirrcili.LetJrcipi=E3+rrci/li·l~Bni2, and (16) can be rewritten as(17)vBrci=Jrcipiq˙pi,where Jrcipi denotes the Jacobian matrix relating the velocity of the upper joint pi to the velocity of the upper leg centroid.In the following section, to obtain the acceleration of the upper leg centroid, (18) can be used, which is given by(18)dl~BniT·lBnidt=-l~BniTl~Bniq˙pilBniT+lBniq˙piTl~BniTl~Bnili.Taking the derivative with respect to time, of (17), and considering (18), the acceleration of the upper leg centroid is given by(19)aBrci=Jrcipi·q¨pi+2l~BniTl~Bniq˙pilBniTq˙pi+lBniq˙piTl~BniTl~Bniq˙pirrcili2.The same procedure may be easily adopted to obtain the velocity and acceleration of the lower leg, which are given by(20)vBtci=rtci·ωai×lBni=Jtcipi·q˙pi,(21)aBtci=Jtcipi·q¨pi-2l~BniTl~Bniq˙pilBniTq˙pi+lBniq˙piTl~BniTl~Bniq˙pirtcili2,where Jtcipi denotes the Jacobian matrix relating the velocity of the upper joint pi to the velocity of the lower leg centroid, Jtcipi=l~BniTl~Bnirtci/li.When the center of mass of the moving platform does not coincide with the origin of the frame{P}, the centrifugal and Coriolis forces will manifest themselves. Therefore the centrifugal and Coriolis terms cannot be ignored when we derive explicit equations for the dynamics of the moving platform. Let pPc denote the position vector of the centroid of the moving platform under body frame {P}, which can be shown by(22)pBc=tB+RPB·pPc.Taking the derivative of both sides of (22) with respect to time, the acceleration of the centroid of the moving platform is given by(23)p¨Bc=t¨B+RPB·p~PcT·RPBT·ω˙+ω~2·RPB·pPc.So far, all the kinematic parameters required for the inverse dynamics are derived. In the next section, attention will be paid to the deduction of dynamic equations of the MMVS.


## 3.1.1. Velocity Analysis

In the base frame, the kinematic equations of theith actuator can be described by(1)lBi=tB+pBi-bBi=tB+RPB·pPi-bBi,where the variable i defines the ith actuator variables, indicating that in general the equations are applicable to any actuator. lBi is the length vector of the actuator with respect to the base frame {B}, tB is the position vector of the body frame, tB=xyzT, and RPB is the rotation matrix of the transformation from the body frame {P} to base frame {B}. The rotation matrix is given by(2)RPB=cos⁡αcos⁡βcos⁡αsin⁡βsin⁡γ-sin⁡αcos⁡γcos⁡αsin⁡βcos⁡γ+sin⁡αsin⁡γsin⁡αcos⁡βsin⁡αsin⁡βsin⁡γ+cos⁡αcos⁡γcos⁡αsin⁡βcos⁡γ-cos⁡αsin⁡γ-sin⁡βcos⁡βsin⁡γcos⁡βcos⁡γ.Taking the derivative of (1) with respect to time, the velocity mapping function can be obtained and is given by(3)q˙pi=l˙Bi=t˙B+ω×RPBpPi=E3RPBp~PiTRPBTq˙,where q˙pi is the velocity of the joint Pi, t˙B and ω are the translational velocity and angular velocity of the moving platform, respectively, t˙B=x˙y˙z˙T, and ω=γ˙β˙α˙T, E3 is a unit 3 × 3 matrix, p~Pi is the skew symmetry matrix of pPi, q˙ is the general velocity of the moving platform, and q˙=t˙BTωTT.LettingJpiq=E3RPBp~PiTRPBT, (3) can be rewritten as(4)q˙pi=Jpiq·q˙,where Jpiq denotes the Jacobian matrix relating the general velocity to the velocity of the upper joint pi.The length of actuator can be described as(5)li=lBiT·lBi.Taking the time derivative of (5), the sliding velocity of the actuator can be obtained:(6)l˙i=lBniT·t˙B+lBniT·ω×RPB·pPi,where lBni is the unit vector of lBi.Equation (6) can be rewritten as(7)l˙i=lBniT,RPB·pPi×lBniT·t˙Bω.Thus (8) may be obtained as(8)L˙=J·q˙,where L˙ denotes actuator sliding velocities, L˙=l˙1,l˙2,l˙3,l˙4,l˙5,l˙6T, and J is the actuator Jacobian matrix expressing mappings from the general velocity to the actuator sliding velocities, which is given by(9)J=lBn1T,RPB·pP1×lBn1TlBn2T,RPB·pP2×lBn2TlBn3T,RPB·pP3×lBn3TlBn4T,RPB·pP4×lBn4TlBn5T,RPB·pP5×lBn5TlBn6T,RPB·pP6×lBn6T.According to their physical meaning, the velocities of the actuator can also be described in terms of the sliding and angular velocities of the actuator underframe{B} as(10)l˙Bi=l˙i·lBni+ωai×lBi,where ωai is the angular velocity of the actuator.Taking the cross product of the above equation withlBni and considering the assumption that no rotation is allowed about the actuator axis (i.e., ωaiT·lBi=0) yield(11)ωai=lBni×l˙Bili=lBni×q˙pili.


## 3.1.2. Acceleration Analysis

For convenience of analysis, the actuator is derived into two parts: the upper leg and the lower leg as shown in Figure5.Figure 5
Diagram of the actuator.Taking the derivative of both sides of (11) with respect to time and taking into account the fact that ωaiT·lBni=0 and εaiT·lBni=0, the angular acceleration of the actuator can be obtained as(12)εai=lBni×q¨pi-2ωai·l˙ili.The translational acceleration of the upper jointPi can be obtained by taking the derivative of (3) with respect to time as described by(13)q¨pi=t¨B+ω˙×RPB·pPi+ω×ω×RPB·pPi.According to their physical meaning, the velocities of the upper jointPi can be described in terms of the velocity of the center of mass of the upper leg of the underframe {B} as(14)q˙pi=vBrci+ωai×rrci·lBni.Rearranging (14) yields(15)vBrci=q˙pi-ωai×rrci·lBni.Substituting (11) into (15), the velocity of the upper leg centroid is given by(16)vBrci=q˙pi+lBni×lBni×q˙pirrcili.LetJrcipi=E3+rrci/li·l~Bni2, and (16) can be rewritten as(17)vBrci=Jrcipiq˙pi,where Jrcipi denotes the Jacobian matrix relating the velocity of the upper joint pi to the velocity of the upper leg centroid.In the following section, to obtain the acceleration of the upper leg centroid, (18) can be used, which is given by(18)dl~BniT·lBnidt=-l~BniTl~Bniq˙pilBniT+lBniq˙piTl~BniTl~Bnili.Taking the derivative with respect to time, of (17), and considering (18), the acceleration of the upper leg centroid is given by(19)aBrci=Jrcipi·q¨pi+2l~BniTl~Bniq˙pilBniTq˙pi+lBniq˙piTl~BniTl~Bniq˙pirrcili2.The same procedure may be easily adopted to obtain the velocity and acceleration of the lower leg, which are given by(20)vBtci=rtci·ωai×lBni=Jtcipi·q˙pi,(21)aBtci=Jtcipi·q¨pi-2l~BniTl~Bniq˙pilBniTq˙pi+lBniq˙piTl~BniTl~Bniq˙pirtcili2,where Jtcipi denotes the Jacobian matrix relating the velocity of the upper joint pi to the velocity of the lower leg centroid, Jtcipi=l~BniTl~Bnirtci/li.When the center of mass of the moving platform does not coincide with the origin of the frame{P}, the centrifugal and Coriolis forces will manifest themselves. Therefore the centrifugal and Coriolis terms cannot be ignored when we derive explicit equations for the dynamics of the moving platform. Let pPc denote the position vector of the centroid of the moving platform under body frame {P}, which can be shown by(22)pBc=tB+RPB·pPc.Taking the derivative of both sides of (22) with respect to time, the acceleration of the centroid of the moving platform is given by(23)p¨Bc=t¨B+RPB·p~PcT·RPBT·ω˙+ω~2·RPB·pPc.So far, all the kinematic parameters required for the inverse dynamics are derived. In the next section, attention will be paid to the deduction of dynamic equations of the MMVS.


## 3.2. Actuator Forces

In this subsection, the Lagrange approach and Newton-Euler method are adopted to develop the dynamic equations of the MMVS. Because of the function of the membrane of the actuators, the effect of gravity on the moving platform and the actuators can be compensated by the elastic forces of the membrane; that is, the moving platform will move to a new equilibrium position; therefore the effect of gravity on the micro-vibration simulator may be ignored.


### 3.2.1. Analysis of a Single Actuator

For the dynamic modeling of the actuator, only the actuator forcesfi and the interaction forces fpi between the upper joints are covered in this subsection. Figure 6 shows a schematic of the actuator. As shown in Figure 6, frame s^i is the body frame of the upper leg which is attached to the center of mass of upper leg with s^1i along theith actuator. Here the over hat indicates a unit of length. The d^i coordinate system is located to the ith actuator at the mass center of the lower leg with d^1i along theith leg.Figure 6
Schematic of the actuator.The Lagrange equation of the actuator can be written as in(24)d∂Ti/∂q˙pidt-∂Ti∂qpi=Qi,where Qi is the generalized force component associated with the generalized coordinate qpi and T is the kinematic energy of the ith actuator, given by(25)Ti=vBrciTmrciBvrci2+vBtciTmtciBvtci2+ωaiTIBrci+IBtciωai2,where IBrci and IBtci are the moments of inertia of the upper leg and lower leg with respect to frame {B}, respectively.The moments of inertiaIBrci and IBtci are discussed in what follows. Assuming that the actuator is axisymmetric and the principal axes of inertia of the upper leg are consistent with the corresponding axes of frame s^i, the principal moments of inertia are denoted by Is^1s^1, Is^2s^2, and Is^3s^3, and the following expression can be obtained:(26)IBrci=RBs^iIs^irciRBs^iT=RBs^i·diag⁡Is^1s^1iIs^2s^2iIs^3s^3i·RBs^iT=Is^1s^1is^1is^1iT+Is^2s^2is^2is^2iT+Is^3s^3is^3is^3iT,where RBs^i is the rotation matrix of the transformation from frame s^i to base frame {B}, given by(27)RBs^i=s^1iT·i^s^2iT·i^s^3iT·i^s^1iT·j^s^2iT·j^s^3iT·j^s^1iT·k^s^2iT·k^s^3iT·k^,where i^=100T, j^=010T, and k^=001T are three unit vectors, which are used to present the base frame {B}.Because the actuator is axisymmetric, the principal moment of inertiaIs^2s^2i is equal to Is^3s^3i. Letting Iri denote Is^2s^2i, (26) can be rewritten by(28)IBrci=Is^1s^1is^1is^1iT+IriE3-s^1is^1iT.Similarly, the moment of inertiaIBtci can be described by(29)IBtci=Id^1d^1id^1id^1iT+ItiE3-d^1id^1iT,where Id^1d^1i is the principal moment of inertia of the lower leg with respect to the axis d^1i of frame d^i and Iti presents the principal moment of inertia of the lower leg with respect to the axis d^2i (or the axis d^3i) of frame d^i.Substituting (17), (20), (28), and (29) into (25) yields(30)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12ωaiTIs^1s^1is^1is^1iT+IriE3-s^1is^1iT+Id^1d^1id^1id^1iT+ItiE3-d^1id^1iTωai.From Figure6, it can be known that s^1i, d^1i, and lni are equal to each other. Thus, (30) can be rewritten as(31)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12Is^1s^1i+Id^1d^1i-Iri-ItiωaiT·lBni2+Iri+ItiωaiT·ωai.Considering the assumption that no rotation is allowed about the actuator axis (i.e.,ωaiT·lBni=0) and substituting (11) into (31), the following expression can be obtained: (32)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12Iri+ItiωaiT·ωai=12q˙piTJrcipimrci+Jtcipimtci+l~BniT·l~BniIri+Itili2q˙pi=12q˙piTMlegiq˙pi,where Mlegi=Jrcipimrci+Jtcipimtci+l~BniT·l~BniIri+Iti/li2.The first and second terms in (24) can be written as(33)d∂T/∂q˙pidt=Mlegiq¨pi+q˙pi2lBniq˙piTl~BniTl~Bnimrcirrcili2+lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT-2lBniTq˙piTl~BniTl~Bni+lBniq˙piTl~BniTl~Bni+l~BniTl~Bniq˙pilBniTmrcirrci2+mtcirtci2+Iri+Itili3,∂T∂qpi=lBniq˙piTq˙pi+2lBniTq˙piq˙pi-3lBniq˙piTlBnilBniTq˙pimrcirrcili2-lBniq˙piTq˙pi+q˙piq˙piTlBni-2lBnilBniTq˙piq˙piTlBnimrcirrci2+mtcirtci2+Iri+Itili3.The right side of (24) can be described as(34)Qi=lBnifi+fpi.Substituting (25), (33), and (34) into (24) yields(35)fpi=Mlegiq¨pi+Caiq˙pi-lBnifi=MlegiJpiqq¨+CaiJpiqq˙+Mlegiω~2RPBpPi-lBnifi,where Cai is a 3×3 matrix, which is given by(36)Cai=mrcirrcili2·lBniq˙piTl~BniTl~Bni+lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT-mrcirrci2+mtcirtci2+Iri+Itili3·lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT.


### 3.2.2. The Completed Dynamic Equations

Assuming that the axial stiffness coefficients of the actuators are equal and lettingk denote the axial stiffness coefficient, then for theith actuator (37)fei=-k·Δli,where fei is the force needed to cause a Δl change in theith actuator length. Assembling the equations for all the actuators and considering (8), (37) becomes(38)fe=-k·J·Δq,where fe is the general elastic force, fe=fe1fe2fe3fe4fe5fe6T.Similarly, the general damping forcefc can be described as (39)fc=-c·J·q˙,where fc is a 6×1 matrix and c is the damping coefficient of the actuator.Based on the virtual work principle, the dynamic equation of the moving platform under the external forces of the loads can be described as(40)ΓT·δq=feT+fcT+FT·δL,where Γ=τ1,τ2,τ3,τ4,τ5,τ6T is the applied generalized force exerted on the moving platform, δq=δx,δy,δz,δγ,δβ,δαT is virtual displacement of the moving platform, F=f1,f2,f3,f4,f5,f6T describes the actuator forces, and δL=δl1,δl2,δl3,δl4,δl5,δl6T gives the virtual displacement of the actuators.(41)δL=J·δq.Substituting (41) into (40), the external force Γ acting on the moving platform can be written as(42)Γ=JTfe+fc+F.Substituting (38) and (39) into (42) yields(43)Γ=-K·Δq-C·q˙+JTF,where K=kJTJ is the matrix for generalized stiffness and C=cJTJ is the matrix for generalized damping.The dynamic model for upper platform as a rigid body ignoring the inertia of the actuators can be derived using the Newton-Euler method [14]:(44)m·p¨Bc=F3×1,RPB·IP·RPBT·ω˙+ω~·RPB·IP·RPBT·ω+m·RPB·pPc×t¨=M3×1,where m is the payload mass, IP is the inertia matrix with respect to frame {P}, F3×1=τ1,τ2,τ3T and M3×1=τ4,τ5,τ6T.Substituting (43) into (44) yields(45)mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBTq¨+000ω~RPB·IP·RPBT+c·JTJq˙+mE30ω~2·RPB·p~Pc+k·JT·J·Δq=JTF.Considering the inertia of the actuators, (45) can be written as(46)mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBTq¨+000ω~RPB·IP·RPBT+c·JTJq˙+mE30ω~2·RPB·p~Pc+k·JT·J·Δq=-HpFP.Expressions forHp and Fp are given by(47)Hp=E3E3E3E3E3E3RPBp~P1RPBTRPBp~P2RPBTRPBp~P3RPBTRPBp~P4RPBTRPBp~P5RPBTRPBp~P6RPBT,FP=fp1Tfp2Tfp3Tfp4Tfp5Tfp6TT.Substituting (35) and (47) into (46) produces the complete dynamic equation of the vibration simulator, which can be described as(48)Mqq¨+Cq,q˙q˙+K·Δq=JTF,where M(q) is a 6×6 mass matrix, Cq,q˙ is a 6×6 matrix of the centrifugal and Coriolis force terms, K is a 6×6 matrix of the generalized stiffness, and F is a 6×1 vector representing actuator forces. They are given by(49)Mq=mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBT+∑i=16JpiqTMlegiJpiq,Cq,q˙q˙=000ω~RPB·IP·RPBT+c·JTJq˙+∑i=16JpiiTCaiJpiiq˙+mE30ω~2RPB·pPc+∑i=16JpiiTMlegiω~2RPB·pPi.Premultiplying both sides of (48) with J-T, the actuator forces are given by(50)F=J-TMqq¨+Cq,q˙q˙+K·Δq.


## 3.2.1. Analysis of a Single Actuator

For the dynamic modeling of the actuator, only the actuator forcesfi and the interaction forces fpi between the upper joints are covered in this subsection. Figure 6 shows a schematic of the actuator. As shown in Figure 6, frame s^i is the body frame of the upper leg which is attached to the center of mass of upper leg with s^1i along theith actuator. Here the over hat indicates a unit of length. The d^i coordinate system is located to the ith actuator at the mass center of the lower leg with d^1i along theith leg.Figure 6
Schematic of the actuator.The Lagrange equation of the actuator can be written as in(24)d∂Ti/∂q˙pidt-∂Ti∂qpi=Qi,where Qi is the generalized force component associated with the generalized coordinate qpi and T is the kinematic energy of the ith actuator, given by(25)Ti=vBrciTmrciBvrci2+vBtciTmtciBvtci2+ωaiTIBrci+IBtciωai2,where IBrci and IBtci are the moments of inertia of the upper leg and lower leg with respect to frame {B}, respectively.The moments of inertiaIBrci and IBtci are discussed in what follows. Assuming that the actuator is axisymmetric and the principal axes of inertia of the upper leg are consistent with the corresponding axes of frame s^i, the principal moments of inertia are denoted by Is^1s^1, Is^2s^2, and Is^3s^3, and the following expression can be obtained:(26)IBrci=RBs^iIs^irciRBs^iT=RBs^i·diag⁡Is^1s^1iIs^2s^2iIs^3s^3i·RBs^iT=Is^1s^1is^1is^1iT+Is^2s^2is^2is^2iT+Is^3s^3is^3is^3iT,where RBs^i is the rotation matrix of the transformation from frame s^i to base frame {B}, given by(27)RBs^i=s^1iT·i^s^2iT·i^s^3iT·i^s^1iT·j^s^2iT·j^s^3iT·j^s^1iT·k^s^2iT·k^s^3iT·k^,where i^=100T, j^=010T, and k^=001T are three unit vectors, which are used to present the base frame {B}.Because the actuator is axisymmetric, the principal moment of inertiaIs^2s^2i is equal to Is^3s^3i. Letting Iri denote Is^2s^2i, (26) can be rewritten by(28)IBrci=Is^1s^1is^1is^1iT+IriE3-s^1is^1iT.Similarly, the moment of inertiaIBtci can be described by(29)IBtci=Id^1d^1id^1id^1iT+ItiE3-d^1id^1iT,where Id^1d^1i is the principal moment of inertia of the lower leg with respect to the axis d^1i of frame d^i and Iti presents the principal moment of inertia of the lower leg with respect to the axis d^2i (or the axis d^3i) of frame d^i.Substituting (17), (20), (28), and (29) into (25) yields(30)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12ωaiTIs^1s^1is^1is^1iT+IriE3-s^1is^1iT+Id^1d^1id^1id^1iT+ItiE3-d^1id^1iTωai.From Figure6, it can be known that s^1i, d^1i, and lni are equal to each other. Thus, (30) can be rewritten as(31)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12Is^1s^1i+Id^1d^1i-Iri-ItiωaiT·lBni2+Iri+ItiωaiT·ωai.Considering the assumption that no rotation is allowed about the actuator axis (i.e.,ωaiT·lBni=0) and substituting (11) into (31), the following expression can be obtained: (32)Ti=12q˙piTJrcipimrci+Jtcipimtciq˙pi+12Iri+ItiωaiT·ωai=12q˙piTJrcipimrci+Jtcipimtci+l~BniT·l~BniIri+Itili2q˙pi=12q˙piTMlegiq˙pi,where Mlegi=Jrcipimrci+Jtcipimtci+l~BniT·l~BniIri+Iti/li2.The first and second terms in (24) can be written as(33)d∂T/∂q˙pidt=Mlegiq¨pi+q˙pi2lBniq˙piTl~BniTl~Bnimrcirrcili2+lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT-2lBniTq˙piTl~BniTl~Bni+lBniq˙piTl~BniTl~Bni+l~BniTl~Bniq˙pilBniTmrcirrci2+mtcirtci2+Iri+Itili3,∂T∂qpi=lBniq˙piTq˙pi+2lBniTq˙piq˙pi-3lBniq˙piTlBnilBniTq˙pimrcirrcili2-lBniq˙piTq˙pi+q˙piq˙piTlBni-2lBnilBniTq˙piq˙piTlBnimrcirrci2+mtcirtci2+Iri+Itili3.The right side of (24) can be described as(34)Qi=lBnifi+fpi.Substituting (25), (33), and (34) into (24) yields(35)fpi=Mlegiq¨pi+Caiq˙pi-lBnifi=MlegiJpiqq¨+CaiJpiqq˙+Mlegiω~2RPBpPi-lBnifi,where Cai is a 3×3 matrix, which is given by(36)Cai=mrcirrcili2·lBniq˙piTl~BniTl~Bni+lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT-mrcirrci2+mtcirtci2+Iri+Itili3·lBniTq˙pil~BniTl~Bni+l~BniTl~Bniq˙pilBniT.


## 3.2.2. The Completed Dynamic Equations

Assuming that the axial stiffness coefficients of the actuators are equal and lettingk denote the axial stiffness coefficient, then for theith actuator (37)fei=-k·Δli,where fei is the force needed to cause a Δl change in theith actuator length. Assembling the equations for all the actuators and considering (8), (37) becomes(38)fe=-k·J·Δq,where fe is the general elastic force, fe=fe1fe2fe3fe4fe5fe6T.Similarly, the general damping forcefc can be described as (39)fc=-c·J·q˙,where fc is a 6×1 matrix and c is the damping coefficient of the actuator.Based on the virtual work principle, the dynamic equation of the moving platform under the external forces of the loads can be described as(40)ΓT·δq=feT+fcT+FT·δL,where Γ=τ1,τ2,τ3,τ4,τ5,τ6T is the applied generalized force exerted on the moving platform, δq=δx,δy,δz,δγ,δβ,δαT is virtual displacement of the moving platform, F=f1,f2,f3,f4,f5,f6T describes the actuator forces, and δL=δl1,δl2,δl3,δl4,δl5,δl6T gives the virtual displacement of the actuators.(41)δL=J·δq.Substituting (41) into (40), the external force Γ acting on the moving platform can be written as(42)Γ=JTfe+fc+F.Substituting (38) and (39) into (42) yields(43)Γ=-K·Δq-C·q˙+JTF,where K=kJTJ is the matrix for generalized stiffness and C=cJTJ is the matrix for generalized damping.The dynamic model for upper platform as a rigid body ignoring the inertia of the actuators can be derived using the Newton-Euler method [14]:(44)m·p¨Bc=F3×1,RPB·IP·RPBT·ω˙+ω~·RPB·IP·RPBT·ω+m·RPB·pPc×t¨=M3×1,where m is the payload mass, IP is the inertia matrix with respect to frame {P}, F3×1=τ1,τ2,τ3T and M3×1=τ4,τ5,τ6T.Substituting (43) into (44) yields(45)mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBTq¨+000ω~RPB·IP·RPBT+c·JTJq˙+mE30ω~2·RPB·p~Pc+k·JT·J·Δq=JTF.Considering the inertia of the actuators, (45) can be written as(46)mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBTq¨+000ω~RPB·IP·RPBT+c·JTJq˙+mE30ω~2·RPB·p~Pc+k·JT·J·Δq=-HpFP.Expressions forHp and Fp are given by(47)Hp=E3E3E3E3E3E3RPBp~P1RPBTRPBp~P2RPBTRPBp~P3RPBTRPBp~P4RPBTRPBp~P5RPBTRPBp~P6RPBT,FP=fp1Tfp2Tfp3Tfp4Tfp5Tfp6TT.Substituting (35) and (47) into (46) produces the complete dynamic equation of the vibration simulator, which can be described as(48)Mqq¨+Cq,q˙q˙+K·Δq=JTF,where M(q) is a 6×6 mass matrix, Cq,q˙ is a 6×6 matrix of the centrifugal and Coriolis force terms, K is a 6×6 matrix of the generalized stiffness, and F is a 6×1 vector representing actuator forces. They are given by(49)Mq=mE3m·RPB·p~PcT·RPBTm·RPB·p~Pc·RPBTRPB·IP·RPBT+∑i=16JpiqTMlegiJpiq,Cq,q˙q˙=000ω~RPB·IP·RPBT+c·JTJq˙+∑i=16JpiiTCaiJpiiq˙+mE30ω~2RPB·pPc+∑i=16JpiiTMlegiω~2RPB·pPi.Premultiplying both sides of (48) with J-T, the actuator forces are given by(50)F=J-TMqq¨+Cq,q˙q˙+K·Δq.


## 3.3. Calculation of Disturbance Forces and Moments

As we know, the object of this study is to develop a micro-vibration simulator which can reproduce desired disturbance forces and moments. So this subsection will calculate the disturbance forces and moments when the actuator forces are applied to the MMVS.The calculation of the disturbance forces and moments is based on the Newton-Euler method. For the convenience of analysis, a frame of reference is attached to the actuator with its origin at the lower joint points,x-axis (hi) along the actuator, y-axis (yi) along the rotation axis (axis fixed to the actuator) of the universal joint, andz-axis (ki) perpendicular to the x- and y-axes according to the right hand rule. Figure 1 shows the frame. Considering the rotational equilibrium of the entire actuator, Euler’s equation for the actuator is obtained as(51)IBrci+IBtciεai+ωai×IBrci+IBtciωai+mrcirrci×aBrci+mtcirtci×aBtci=lBi×fpi+Muihi,where rrci=(li-rrci)lBni, rtci=rtcilBni, and Mui is the magnitude of the constraint moment at the universal joint acting about the hi axis.Taking cross products of both sides of (51) with hi, the scalar unknown Mui is given by (52)Mui=lBniTCilBniThi,where Ci is a 3×1 matrix, which is given by(53)Ci=IBrci+IBtciεai+ωai×IBrci+IBtciωai+mrcirrci×aBrci+mtcirtci×aBtci.Ignoring the damping forces of the actuator and considering the upper leg of the actuator, Newton’s equation can be given by(54)fi+lBniTfpi-k·Δli-mrcilBniT·aBrci=0.Solving (54), the constraint force at the spherical joint can be rewritten as(55)fpi=XilBni+Di.ExpressionsXi and Di are given by(56)Xi=mrcilBniT·aBrci+k·Δli-fi,Di=Ci×lBni-Muihi×lBnili.Considering the translational equilibrium, Newton’s equation for the entire actuator can be written as(57)mrciBarci+mtciBatci-fpi=fui,where fui is the constraint force at the universal joint.Substituting (55) into (57), the constraint force fui is given by(58)fui=mrciBarci+mtciBatci-XilBni+Di.Consequently, the disturbance forces acting on the center of the bottom of the base plane, which are produced by the simulator, can be written as(59)Fm=∑i=16fui.Similarly, the disturbance moments acting on the center bottom of the base plane, produced by the simulator, which are shown in Figure6, can be written as(60)Mm=∑i=16rOBi×fui+Muihi,where rOBi is the vector from the origin Og to the lower joint point Bi.


## 3.4. Planning of Disturbance Forces and Moments

In the last subsection, the relationships between the disturbance forces and moments and the actuator forces are derived. This subsection will deduce the relationship between the deduced forces and moments and the generalized acceleration of the moving platform and then plan the disturbance forces and moments.Because micro-vibrations have the characteristics of a wide frequency range and small amplitudes, a few approximations can be adopted to plan the disturbance forces and moments. For example, the third term on the right side of (13) is a square term of angular velocity, so it can be ignored when the angular velocity of the moving platform is small. The simplified acceleration of the upper jointPi can be described as(61)p¨Bi≅t¨B+ω˙×RPB·pPi≅E3p~BiTt¨Bω˙=Tpiq¨,where Tpi=E3p~BiT.Similarly, ignoring the second terms of (19) and (21), respectively, the acceleration of the upper and lower legs is given by (62)aBrci=Jrcipi·q¨pi,aBtci=Jtcipi·q¨pi.Substituting (61) into (62) yields(63)aBrci=JrcipiE3p~BiTt¨Bω˙=JrcipiTpiq¨,aBtci=JtcipiE3p~BiTt¨Bω˙=JtcipiTpiq¨.The angular velocity of the actuator is very small, so it can also be ignored. The angular acceleration of the actuator can be described as(64)εai≅l~BniaBpili.After obtaining the simplified expressions presented as (61), (63), and (64), the following expressions are obtained:(65)Xi≅mrcilBniT·JrcipiTpiq¨+k·Δli-fi,Ci≅IBrci+IBtci·l~Bnili+mrcir~rciJrcipi+mtcir~tciJtcipiTpiq¨,Di≅-l~BniliIBrci+IBtci/li·l~Bni+mrcir~rciJrcipi+mtcir~tciJtcipiTpiq¨=Λiq¨,where Λi=-l~Bni/liIBrci+IBtci/li·l~Bni+mrcir~rciJrcipi+mtcir~tciJtcipiTpi.Therefore, the constraint forces at the spherical and universal joints can be rewritten as(66)fpi=XilBni+Λiq¨,fui=mrciJrcipi+mtciJtcipi-mrcilBnilBniTJrcipiTpi-Λi+lBnieiJ-TMqq¨=Θuiq¨,where ei=[0⋯1⋯0] is a 1×6 matrix, theith element is 1, and the other elements are 0, and Θui is 3×1 matrix, Θui=mrciJrcipi+mtciJtcipi-mrcilBnilBniTJrcipiTpi-Λi+lBnieiJ-TMq.The disturbance forces and moments acting on the center of the base plane’s bottom can be rewritten as(67)Fm=∑i=16fui=∑i=16Θuiq¨=ΘTq¨,Mm=∑i=16rOBi×fui+Muihi=∑i=16r~OBiΘuiq¨=ΘRq¨.Equation (67) can be written in a matrix form as(68)FmMm=ΘTΘRq¨.Taking both sides of (68) with ΘTTΘRT-T, the generalized acceleration of the moving platform can be given by(69)q¨=ΘTΘR-1FmMm.LetFmd and Mmd denote the desired disturbance forces and moments, respectively; the desired generalized acceleration of the moving platform can be described by(70)q¨d=ΘTΘR-1FmdMmd.Assuming that the generalized velocities and generalized displacements of the moving platform are, respectively, the time integral and time double integral of the generalized acceleration, the following equations can be obtained:(71)q˙d=∫q¨ddt,qd=∬q¨ddt.Substituting (70) and (71) into (50) yields the desired actuator forces, which are given by(72)Fd=J-TMqdq¨d+Cqd,q˙dq˙d+Kqdqd.So far, the relations between the actuator forces and the desired disturbance forces and moments are obtained.


## 4. Numerical Simulation and Analysis

### 4.1. Normal Mode Analysis

It is very important to consider the dynamic characteristics at an earlier design stage, such as the natural frequencies. As described in (48), the whole dynamic equations of the MMVS are developed. The Coriolis and centrifugal force terms as the square terms of velocity may be neglected when the velocity is small. Thus, the free vibration of the MMVS may be described as(73)Mqq¨+Cq,q˙q˙+Kqq=0.It should be noted that the mass matrixM(q) is given in a concise matrix form, which considers the impact of the payload decent ration and mass of the actuator. The natural frequencies and mode shapes (or eigenvectors) of the MMVS can be obtained by solving this equation with numerical methods. The structural parameters and mass properties of the MMVS are shown in Tables 1 and 2, respectively.Table 1
Structural parameters of the multi-degree-of-freedom micro-vibration simulator.
Notation
Specification
Value
RP
Moving platform radius
0.127 m
Rb
Base plane radius
0.18 m
H
Height of the origin of the body frame in the base frame
0.22 m
φ
Moving platform central angle
30°
θ
Base plane central angle
90°
rrci
Length from the upper joint to the centroid of the upper leg of the actuator
85.9 mm
rtci
Length from the lower joint to the centroid of the lower leg of the actuator
84.8 mm
zcm
Height of the centroid of the moving platform in the body frame
27.8 mm
k
Axial stiffness of membrane
1.6 × 104 N m−1Table 2
Mass properties of the multi-degree-of-freedom micro-vibration simulator.
Notation
Specification
Value
m
Mass of the moving platform
4.934 kg
Ixx
Moment of inertia of the moving platform aboutx-axis
3.0 × 10−2 kg m2
Iyy
Moment of inertia of the moving platform abouty-axis
3.0 × 10−2 kg m2
Izz
Moment of inertia of the moving platform aboutz-axis
5.9 × 10−2 kg m2
mrci
Mass of the upper leg of the actuator
0.443 kg
mtci
Mass of the lower leg of the actuator
1.797 kg
Iri
Moment of inertia of the upper leg
9.0 × 10−4 kg m2
Iti
Moment of inertia of the lower leg
2.3 × 10−3 kg m2The finite element method is also adopted to calculate the natural frequencies and mode shapes of the MVVS. The aim here is to compare and cross-check the theoretically obtained results with a totally different method (in this case a NASTRAN® simulation). The reference MVVS modeled in the simulation has exactly the same characteristics as it does in the theoretical model. Table3 shows a comparison between the six natural frequencies obtained by the analytical equations formulated in this paper and the simulation. The comparison shows an excellent agreement between the results in terms of modal properties.Table 3
Natural frequencies obtained by the developed analytical model and the simulation.
Method
f
1
f
2
f
3
f
4
f
5
f
6
(Hz)
(Hz)
(Hz)
(Hz)
(Hz)
(Hz)
Theoretical
4.79
4.79
7.33
16.43
17.45
17.45
Simulation
4.86
4.86
7.32
15.99
16.73
16.73
Relative error
1.5%
1.5%
0.1%
2.7%
4.3%
4.3%The mode shapes in the work space directly represent the position and orientation of the moving platform corresponding to each natural frequency. Substituting the design parameters and natural frequencies of the MVVS into (73), the modal matrix is obtained using MATLAB and rearranged according to the natural frequencies sequence as described with expression (74). Having normalized the modal matrix, the mode shapes obtained by the theoretical method and the simulation are exactly in correspondence; see Figure 7.(74)ϕ=10000-0.010-100-0.0100001000-0.8900-100.890000100.021000.Figure 7
The first six mode shapes of the simulator: (a) shearing mode alongx at 4.86 Hz, (b) shearing mode along y at 4.86 Hz, (c) twisting mode around z at 7.32 Hz, (d) piston mode along z at 15.99, (e) tilting mode around x at 16.73 Hz, and (f) tilting mode around y at 16.73 Hz.
(a)
(b)
(c)
(d)
(e)
(f)


### 4.2. Cosimulation Verification

In this section, a cosimulation using ADAMS and MATLAB/Simulink is adopted to verify the validity of the dynamic model and the feasibility of the MMVS [22, 23]. The control system is constructed in the MATLAB/Simulink environment and the virtual prototype of the MMVS is built in ADAMS. The six actuator forces of the actuators are used as inputs to drive the ADAMS model and the disturbance forces and moments of the base plane acting on the ground are used as the outputs of the ADAMS plant model. Figure 8 shows the simplified ADAMS model and the cosimulation block diagram. The corresponding Simulink model is shown in Figure 9. The model has seven blocks on the top level. The first, leftmost block is the command generator block which can generate the desired disturbance forces and moments that we design. The second is the configuration parameters block used to assign the structural parameters of the MMVS. The third block is used to assign the actuator parameters, and the third block is used to assign the stiffness coefficient of the membrane, damping coefficient of the single actuator, and the mass properties of payload. The desired actuator forces are obtained from the output of the fifth block which is the MATLABFcn block including the dynamic equations of the MMVS. The sixth is a subsystem block representing the virtual prototype building in ADAMS. The seventh is used to plot the actual disturbance forces and moments measured by “sensors” in ADAMS.Figure 8
Cosimulation block diagram.Figure 9
Simulink model of cosimulation of the MMVS.According to the initial design goal, the MMVS needs to generate the frequency spectrum of the disturbance forces and moments produced by real flight R/MWAs. Table4 shows the target disturbance forces and moments, which are the experimental data reported in [4].Table 4
Target disturbance forces and moments.
Frequency (Hz)
F
x (N)
F
y (N)
F
z (N)
M
x (N·m)
M
y (N·m)
M
z (N·m)
50
0.1
0.1
0
0.05
0.05
0
100
0.2
0.5
0
0.15
0.15
0
200
0.2
0.2
1.2
0.05
0.05
0
250
1.8
1.8
3.2
0.4
0.38
0.05Figures10 and 11 show the disturbance forces and moments produced by the MMVS in the time and frequency domains, respectively. From Figure 11, it is clear that the disturbance forces and moments generated by the simulator are quite consistent with the target MWA disturbance data reported in [4]. Moreover, there is no coupling between the generated disturbance forces and moments. This indicates that the micro-vibration simulator can reproduce the disturbance forces and moments along any arbitrary direction. Table 5 shows the frequency spectrum difference between the MWA disturbance and the simulator disturbance. As shown in Table 4, the frequency components of their disturbance forces and moments are identical. The greatest differences in the disturbance force amplitude and the disturbance moment amplitude are 4.19% and 7.5%, respectively. The smallest differences in the disturbance force amplitude and the disturbance moment amplitude are 0.08% and 1.68%, respectively.Table 5
Frequency spectrum difference between MWA disturbance and simulator disturbance.
Frequency (Hz)
Difference (%)
F
x
F
y
F
z
M
x
M
y
M
z
50 Hz
3.20
4.19
0.00
2.56
4.42
0.00
100 Hz
1.25
0.08
0.00
3.67
2.73
0.00
200 Hz
0.35
0.40
0.58
7.30
2.42
0.00
250 Hz
0.44
0.50
0.44
7.50
3.29
1.68Figure 10
Disturbance forces and moments in the time domain: (a) disturbance forces along thex-axis, (b) disturbance forces along the y-axis, (c) disturbance forces along the z-axis, (d) disturbance moments around the x-axis, (e) disturbance moments around the y-axis, and (f) disturbance moments around the z-axis.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Disturbance forces and moments in the frequency domain: (a) disturbance forces along thex-axis, (b) disturbance forces along the y-axis, (c) disturbance forces along the z-axis, (d) disturbance moments around the x-axis, (e) disturbance moments around the y-axis, and (f) disturbance moments around the z-axis.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.1. Normal Mode Analysis

It is very important to consider the dynamic characteristics at an earlier design stage, such as the natural frequencies. As described in (48), the whole dynamic equations of the MMVS are developed. The Coriolis and centrifugal force terms as the square terms of velocity may be neglected when the velocity is small. Thus, the free vibration of the MMVS may be described as(73)Mqq¨+Cq,q˙q˙+Kqq=0.It should be noted that the mass matrixM(q) is given in a concise matrix form, which considers the impact of the payload decent ration and mass of the actuator. The natural frequencies and mode shapes (or eigenvectors) of the MMVS can be obtained by solving this equation with numerical methods. The structural parameters and mass properties of the MMVS are shown in Tables 1 and 2, respectively.Table 1
Structural parameters of the multi-degree-of-freedom micro-vibration simulator.
Notation
Specification
Value
RP
Moving platform radius
0.127 m
Rb
Base plane radius
0.18 m
H
Height of the origin of the body frame in the base frame
0.22 m
φ
Moving platform central angle
30°
θ
Base plane central angle
90°
rrci
Length from the upper joint to the centroid of the upper leg of the actuator
85.9 mm
rtci
Length from the lower joint to the centroid of the lower leg of the actuator
84.8 mm
zcm
Height of the centroid of the moving platform in the body frame
27.8 mm
k
Axial stiffness of membrane
1.6 × 104 N m−1Table 2
Mass properties of the multi-degree-of-freedom micro-vibration simulator.
Notation
Specification
Value
m
Mass of the moving platform
4.934 kg
Ixx
Moment of inertia of the moving platform aboutx-axis
3.0 × 10−2 kg m2
Iyy
Moment of inertia of the moving platform abouty-axis
3.0 × 10−2 kg m2
Izz
Moment of inertia of the moving platform aboutz-axis
5.9 × 10−2 kg m2
mrci
Mass of the upper leg of the actuator
0.443 kg
mtci
Mass of the lower leg of the actuator
1.797 kg
Iri
Moment of inertia of the upper leg
9.0 × 10−4 kg m2
Iti
Moment of inertia of the lower leg
2.3 × 10−3 kg m2The finite element method is also adopted to calculate the natural frequencies and mode shapes of the MVVS. The aim here is to compare and cross-check the theoretically obtained results with a totally different method (in this case a NASTRAN® simulation). The reference MVVS modeled in the simulation has exactly the same characteristics as it does in the theoretical model. Table3 shows a comparison between the six natural frequencies obtained by the analytical equations formulated in this paper and the simulation. The comparison shows an excellent agreement between the results in terms of modal properties.Table 3
Natural frequencies obtained by the developed analytical model and the simulation.
Method
f
1
f
2
f
3
f
4
f
5
f
6
(Hz)
(Hz)
(Hz)
(Hz)
(Hz)
(Hz)
Theoretical
4.79
4.79
7.33
16.43
17.45
17.45
Simulation
4.86
4.86
7.32
15.99
16.73
16.73
Relative error
1.5%
1.5%
0.1%
2.7%
4.3%
4.3%The mode shapes in the work space directly represent the position and orientation of the moving platform corresponding to each natural frequency. Substituting the design parameters and natural frequencies of the MVVS into (73), the modal matrix is obtained using MATLAB and rearranged according to the natural frequencies sequence as described with expression (74). Having normalized the modal matrix, the mode shapes obtained by the theoretical method and the simulation are exactly in correspondence; see Figure 7.(74)ϕ=10000-0.010-100-0.0100001000-0.8900-100.890000100.021000.Figure 7
The first six mode shapes of the simulator: (a) shearing mode alongx at 4.86 Hz, (b) shearing mode along y at 4.86 Hz, (c) twisting mode around z at 7.32 Hz, (d) piston mode along z at 15.99, (e) tilting mode around x at 16.73 Hz, and (f) tilting mode around y at 16.73 Hz.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.2. Cosimulation Verification

In this section, a cosimulation using ADAMS and MATLAB/Simulink is adopted to verify the validity of the dynamic model and the feasibility of the MMVS [22, 23]. The control system is constructed in the MATLAB/Simulink environment and the virtual prototype of the MMVS is built in ADAMS. The six actuator forces of the actuators are used as inputs to drive the ADAMS model and the disturbance forces and moments of the base plane acting on the ground are used as the outputs of the ADAMS plant model. Figure 8 shows the simplified ADAMS model and the cosimulation block diagram. The corresponding Simulink model is shown in Figure 9. The model has seven blocks on the top level. The first, leftmost block is the command generator block which can generate the desired disturbance forces and moments that we design. The second is the configuration parameters block used to assign the structural parameters of the MMVS. The third block is used to assign the actuator parameters, and the third block is used to assign the stiffness coefficient of the membrane, damping coefficient of the single actuator, and the mass properties of payload. The desired actuator forces are obtained from the output of the fifth block which is the MATLABFcn block including the dynamic equations of the MMVS. The sixth is a subsystem block representing the virtual prototype building in ADAMS. The seventh is used to plot the actual disturbance forces and moments measured by “sensors” in ADAMS.Figure 8
Cosimulation block diagram.Figure 9
Simulink model of cosimulation of the MMVS.According to the initial design goal, the MMVS needs to generate the frequency spectrum of the disturbance forces and moments produced by real flight R/MWAs. Table4 shows the target disturbance forces and moments, which are the experimental data reported in [4].Table 4
Target disturbance forces and moments.
Frequency (Hz)
F
x (N)
F
y (N)
F
z (N)
M
x (N·m)
M
y (N·m)
M
z (N·m)
50
0.1
0.1
0
0.05
0.05
0
100
0.2
0.5
0
0.15
0.15
0
200
0.2
0.2
1.2
0.05
0.05
0
250
1.8
1.8
3.2
0.4
0.38
0.05Figures10 and 11 show the disturbance forces and moments produced by the MMVS in the time and frequency domains, respectively. From Figure 11, it is clear that the disturbance forces and moments generated by the simulator are quite consistent with the target MWA disturbance data reported in [4]. Moreover, there is no coupling between the generated disturbance forces and moments. This indicates that the micro-vibration simulator can reproduce the disturbance forces and moments along any arbitrary direction. Table 5 shows the frequency spectrum difference between the MWA disturbance and the simulator disturbance. As shown in Table 4, the frequency components of their disturbance forces and moments are identical. The greatest differences in the disturbance force amplitude and the disturbance moment amplitude are 4.19% and 7.5%, respectively. The smallest differences in the disturbance force amplitude and the disturbance moment amplitude are 0.08% and 1.68%, respectively.Table 5
Frequency spectrum difference between MWA disturbance and simulator disturbance.
Frequency (Hz)
Difference (%)
F
x
F
y
F
z
M
x
M
y
M
z
50 Hz
3.20
4.19
0.00
2.56
4.42
0.00
100 Hz
1.25
0.08
0.00
3.67
2.73
0.00
200 Hz
0.35
0.40
0.58
7.30
2.42
0.00
250 Hz
0.44
0.50
0.44
7.50
3.29
1.68Figure 10
Disturbance forces and moments in the time domain: (a) disturbance forces along thex-axis, (b) disturbance forces along the y-axis, (c) disturbance forces along the z-axis, (d) disturbance moments around the x-axis, (e) disturbance moments around the y-axis, and (f) disturbance moments around the z-axis.
(a)
(b)
(c)
(d)
(e)
(f)Figure 11
Disturbance forces and moments in the frequency domain: (a) disturbance forces along thex-axis, (b) disturbance forces along the y-axis, (c) disturbance forces along the z-axis, (d) disturbance moments around the x-axis, (e) disturbance moments around the y-axis, and (f) disturbance moments around the z-axis.
(a)
(b)
(c)
(d)
(e)
(f)


## 5. Conclusions

In this study, a multi-degree-of-freedom micro-vibration simulator (MMVS) based on the Gough-Stewart platform, which can generate the disturbance frequency spectrum of a real flight reaction/momentum wheel assemblies (R/MWAs), is proposed. The complete inverse kinematic and dynamic equations of the simulator were derived using the Newton-Euler method and Lagrange approach. These kinematic and dynamic models were proved to be valid. By approximation of a feasible model, a disturbance planning method to calculate the excitation forces of the actuators was deduced. To verify the feasibility of this micro-vibration simulator, cosimulation was adopted. The simulation results show that the simulator induced disturbances were almost identical to those of a real flight R/MWAs, and there was no coupling between the disturbance forces and moments.In this study, valuable findings for a micro-vibration simulator of the R/MWAs were obtained. However, the dynamic models of the simulator were derived based on considering the simulator as a rigid system. In future studies, the structure resonance caused by flexing will be taken into account while planning the excitation forces of the actuators.



---

*Source: 4840514-2017-08-02.xml*