# Research on Microvibrations Generated by a Control Moment Gyroscope on a Flexible Interface Based on a Dynamic Substructure Method

Authors: Xiongfei Li; Wei Cheng
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

Microvibrations generated by a control moment gyroscope (CMG) will couple with a flexible spacecraft structure, and this seriously degrades certain point performances of a spacecraft. This study focuses on investigating the coupled microvibrations caused by a CMG on a flexible interface by using a dynamic substructure method (DSM). First, a DSM based on a frequency response function (FRF) is established, and this method is used to simultaneously synthesize multiple substructures with different coordinates irrespective of whether their connection interfaces are rigid or flexible. Second, the bearings and the gimbal servo system are simplified into linear springs, and therefore the CMG model is equivalent to a mass-spring-damping system with eighteen degrees of freedom (DOFs). Third, the established DSM is employed to deduce the FRF matrix of the CMG-flexible interface coupling system that consists of CMG mounted on an aluminum honeycomb sandwich palate (AHSP). Dynamic responses of the coupling system are calculated by the derived FRF matrix. Finally, MATLAB and multibody dynamics simulations are conducted to analyze and validate the dynamic responses of the coupling system obtained by the DSM. The results indicate that the DSM is appropriate to predict the coupled microvibrations of CMG on a flexible interface and exhibits high prediction accuracy and computational efficiency.


---

## Body

## 1. Introduction

With the development of high-precision observation spacecrafts, spacecraft microvibrations are increasingly important in the design process of spacecrafts and result in stringent requirements for the spacecraft platform stability [1, 2]. Microvibrations of spacecrafts are characterized by a low-amplitude and wide-frequency band, and they significantly influence both imaging quality and pointing accuracy of a spacecraft. They are typically caused by internal mechanisms on board spacecrafts, such as the control moment gyroscope (CMG), reaction wheel assembly (RWA), cryocooler assembly (CCA), camera shutter assembly (CSA), and solar array drive assembly (SADA), which in this situation are termed as microvibration sources [3]. Among these mechanisms, RWA and CMG are always considered the most prominent disturbance sources. Their disturbances are distributed in middle-to-high-frequency band and are difficult to control and reduce [4, 5].Additionally, RWA and CMG are commonly used as actuators of attitude control in spacecrafts via the generation of reaction torque by high-speed rotary flywheels although their operating principles are different. RWA provides control torque for spacecrafts via adjusting the rotational acceleration of the flywheel while CMG generates gyroscopic torque by turning the angular momentum vector of the flywheel in which the rotating speed is constant and considerably higher than the maximum of RWA and is as high as approximately 6000–10,000 rpm [6]. The rotary flywheel also provides adverse vibrations while providing control torque. These vibrations are modulated and magnified by internal dynamic mechanisms and transmitted through the spacecraft structure, exciting modes of the spacecraft structure, and significantly affect its performance [7]. The flywheel disturbances are mainly due to the following four factors: flywheel mass imbalance, internal resonance, imperfection in mechanical bearings, and motor ripple, in which flywheel mass imbalance is almost considered the major factor.Several studies investigated the disturbances induced by the rotary flywheel [8–17]. Most studies focus on disturbances produced by RWA. Reference [8] proposed an empirical model of RWA that assumes that RWA disturbances consist of a series of harmonics at discrete frequencies with amplitude that is proportional to the square of the flywheel speed. However, the empirical model fails to express the amplification when the disturbance harmonics cross structural modes. Therefore, [9] developed an analytical model that captures the dynamic interaction between the structural modes and the inertia properties of the flywheel and gyroscopic stiffening effects. The aforementioned disturbance models of RWA are based on the rigid interface condition assumption. However, in terms of application, RWA is mounted on a flexible spacecraft interface, and it is highly likely that dynamic coupling between the RWA and spacecraft interface is induced during operation. The RWA disturbances induce vibrations in the spacecraft and excite its flexible modes. The spacecraft vibrations subsequently drive the RWA and create additional disturbances, and these increase the complexity of RWA microvibrations [9]. References [10–13] investigated the coupled disturbances between the RWA and flexible interface by using dynamic mass measurement techniques. This method uses a “force filter” to correct the RWA disturbance test spectrums of the rigid interface and applies the corrected spectra to a spacecraft frequency response function (FRF) to predict its performance. References [14–17] studied coupled microvibrations of a cantilever-configured RWA with a flexible platform based on theoretical modeling and experimental measurement. A seismic microvibration measurement system (SMVMS) was developed to measure the coupled microvibrations, and an analytical model of the RWA and the test system was established to analyze the coupled microvibrations in which the importance of the RWA driving point accelerance was thoroughly discussed.When compared to RWA, studies for CMG disturbances are relatively less. References [18, 19] established the disturbance models of CMG on a rigid interface via the Lagrange energy method and Euler equation, respectively. Both models were verified by experimental tests. In a manner similar to RWA, the microvibration characteristics of CMG on the flexible interface are more complicated. Reference [20] analyzed the influence of installation bracket stiffness on CMG microvibration properties by testing CMG microvibrations on two installation brackets with different stiffness values. Reference [21] investigated the effect of the flexible interface on CMG microvibrations by comparing microvibration testing results of CMG on rigid and flexible interfaces. However, the aforementioned studies did not deeply discuss the coupled microvibration characteristics of CMG on the flexible interface because only a few effective methods exist to describe and analyze the coupled microvibrations.Because of the structural complexity of the coupled system, this study adopts a dynamic substructure method (DSM) to investigate the coupled microvibration characteristics. The basic idea of DSM involves decomposing an assembly system into several components and predicting the dynamic properties of the assembly structure from its components. The DSM is divided into two groups, namely, one in which modal data is used and another in which FRF is used directly [22, 23]. The FRF-using DSM is selected in this study. The first FRF-using DSM is the impedance coupling method (IC method) that uses the receptance matrix of substructures to obtain that of the assembly structure [24]. However, this method is computationally inefficient because it should initially inverse the FRF matrix of each substructure. Reference [25] proposed a receptance coupling method (RC method) that is more computationally efficient by synthesizing the FRF of substructures directly although it is only used for coupling two independent substructures. Reference [26] presented a generalized receptance coupling method (GRC method) that is used to couple several substructures together in a single step. Reference [27] developed a general joint description method (GJDM) that is very similar to the GRC method in terms of formulation and convenience to resolve the effect of elastic media between substructures. In this study, a DSM convenient for synthesizing multiple different coordinate substructures with hybrid connections is developed. Subsequently, this method is employed to deduce the FRF matrix of the CMG-flexible interface coupling system that is used to calculate the dynamic responses of the coupling system. Finally, MATLAB simulations are conducted to analyze the dynamic responses of the coupling system, and multibody dynamics simulations are performed to analyze and verify the accuracy of the coupling disturbance analysis method based on the DSM.


## 2. Dynamic Substructure Method

### 2.1. Synthesis Criterions of the Flexible Connection

Figure1 describes a substructure system divided from an assembly substructure with intricate connections that contains both flexible and rigid connections. Each connection of the substructure system is considered a substructure and is defined as a connection substructure (CS). A rigid connection essentially belongs to a special flexible connection in which stiffness tends towards infinity. Therefore, all the connections are considered flexible connections in the following derivation.Figure 1
Substructure system with intricate connections.We define two groups of connection coordinates in which a group corresponds to connection coordinatesb′, b¯, and b located on connection substructures, substructures, and assembly structure, respectively, and they are on a side of all the connections corresponding to another group, namely, c′, c¯, and c on another side. For a multidegree of freedom (DOF) connection substructure j, each DOF can be equivalent to a spring with equivalent stiffness-damping, as shown in Figure 2, in which each full line in the CS represents a DOF and the dotted line represents the remaining DOFs. The relation between response and excitation of the jth connection substructure is determined as follows:
(1)Fbj′Fcj′=Zj−Zj−ZjZjXbj′Xcj′,Zj=Kj+sCj,where j=1, 2,…,n, n is the sum of the connection substructures, Zj=zabo×o refers to the connection impedance matrix of the jth connection substructure, o denotes the DOF sum of the jth connection substructure, and zab denotes the connection impedance between ath and bth DOFs that denotes the force needed to act on the bth DOF when the ath DOF happens to be a unit displacement. When no couple exists between different connection DOFs, Zj corresponds to a diagonal matrix. Additionally, Kj and Cj represent the equivalent connection stiffness and damping matrixes, respectively, of the jth connection substructure; s=wi=2πfi, where i denotes the square root of −1 and w and f denote the vibration angular frequency and vibration frequency, respectively.Figure 2
Connection substructure between two substructures.Furthermore,Rj is defined as the connection receptance matrix of the jth connection substructure that is the inverse matrix of the connection impedance matrix and is calculated by (2). For a rigid connection, the connection receptance matrix is considered a zero matrix.
(2)Rj=Zj−1=Kj+sCj−1.Subsequently, the relation between response and excitation of thejth connection substructure is expressed by the connection receptance matrix as follows:
(3)Xbj′−Xcj′=RjFbj′=−RjFcj′.Similarly, the relation between response and excitation of all the connection substructures is expressed as follows:(4)Xb′−Xc′=RFb′=−RFc′,where R denotes the connection receptance matrix of the substructure system, and it is obtained by connection receptance matrixes of all the connection substructures in the diagonal form and expressed as follows:
(5)R=diagR1R2⋯Rn.


### 2.2. Formula Derivation

For the substructure system shown in Figure1, FRF matrixes of all the substructures are partitioned in internal coordinates and connection coordinates b¯ and c¯ that denote the relationship between response and excitation and are expressed as follows:
(6)Xa¯Xb¯Xc¯=H¯Fa¯Fb¯Fc¯=Ha¯a¯Ha¯b¯Ha¯c¯Hb¯a¯Hb¯b¯Hb¯c¯Hc¯a¯Hc¯b¯Hc¯c¯Fa¯Fb¯Fc¯.Accordingly, the relationship between response and excitation of the assembly structure is expressed as follows:(7)XaXbXc=HFaFbFc=HaaHabHacHbaHbbHbcHcaHcbHccFaFbFc,where subscripts a¯ and a represent the internal coordinates of the substructure system and the assembly structure, respectively.The basic conditions in the process of substructure synthesis are compatibility and equilibrium that are given as follows:(8)Xb¯Xc¯=Xb′Xc′=XbXc,(9)Fb¯Fc¯+Fb′Fc′=FbFc.In the substructure synthesis process, the response and excitation of the internal coordinates are not changed, and they are expressed as follows:(10)Xa¯=Xa,(11)Fa¯=Fa.By combining (4), (6), and (8), we obtain the following expression:
(12)RFb′=−RFc′=Hb¯a¯−Hc¯a¯Hb¯b¯−Hc¯b¯Hb¯c¯−Hc¯c¯TFa¯Fb¯Fc¯.The expression forFb′ and Fc′ is expressed as (13) by substituting (9) and (11) into (12) as follows:
(13)Fb′=−Fc′=R+Hb¯b¯+Hc¯c¯−Hc¯b¯−Hb¯c¯−1Hb¯a¯−Hc¯a¯Hb¯b¯−Hc¯b¯Hb¯c¯−Hc¯c¯TFaFbFc.Using (4), (8), (9), and (11), the relationship between response and excitation of the substructure system given by (6) is reexpressed as follows:
(14)XaXbXc=H¯FaFbFc+0−Fb′Fb′.By substituting (13) into (14) for Fb′, the FRF matrix of the assembly structure is obtained and expressed as (15) in which the FRF matrix of any assembly structure is calculated by the FRF and connection receptance matrixes of its substructures.
(15)H=H¯−Ha¯b¯−Ha¯c¯Hb¯b¯−Hb¯c¯Hb¯c¯−Hc¯c¯R+Hb¯b¯+Hc¯c¯−Hc¯b¯−Hb¯c¯−1Hb¯a¯−Hc¯a¯Hb¯b¯−Hc¯b¯Hb¯c¯−Hc¯c¯T.


### 2.3. Process of Substructure Synthesis

We consider a three-substructure system an example, and the process of substructure synthesis is demonstrated in this section. Figure3 shows the three-substructure system where substructure II connects to substructures I and III via connection substructures 1 and 2, respectively. Coordinates o1x1y1z1 (coordinate 1), o2x2y2z2 (coordinate 2), and o3x3y3z3 (coordinate 3) are the local coordinates of the three substructures, respectively, and are discordant from each other.Figure 3
Substructure system with three substructures.The FRF matrixes of the three substructures in their local coordinates are defined as follows:(16)HI=HaaIHab1IHb1aIHb1b1I,HII=HaaIIHac1IIHac2IIHc1aIIHc1c1IIHc1c2IIHc2aIIHc2c1IIHc2c2II,HIII=HaaIIIHab2IIIHb2aIIIHb2b2III.The FRF matrix of the substructure system is obtained by distributing FRF matrixes of the three substructures in internal coordinatea and connection coordinates b and c and expressed as follows:
(17)Due to the inconformity of the three coordinates, the connection coordinatesb and c are not in agreement. The FRF matrix in (17) is substituted into (15) only if the connection coordinates b and c are transformed into accordant coordinates. It is assumed that the connection coordinates of substructures I and III are accordant to those of substructure II after the transformation, and thus it is necessary to transform the connection coordinates of substructures I and III by a connection coordinate transformation matrix that is defined as follows:
(18)A=diagIuA1A2Iv,where I denotes the unit matrix, subscripts u and v refer to their order that is the sum of DOFs of all the internal and connection coordinates, and matrixes A1 and A2 represent the orientation cosine matrixes of all the connection coordinate DOFs of substructure II with those of substructures I and III, respectively.After the transformation via the connection coordinate transformation matrix, the FRF matrix of the substructure system is reexpressed as follows:(19)The connection receptance matrixes of connection substructures 1 and 2 are defined asR1 and R2, respectively, in coordinate 2. With respect to (5), the connection receptance matrix of the substructure system is expressed as follows:
(20)R=R100R2.By substituting (19) and (20) into (15), the FRF matrix of the assembly structure is achieved. The steps involved in the substructure synthesis process can be summarized as follows:
(1)
Dividing a complex assembly structure into several substructures depending on the analysis demand
(2)
Calculating the FRF matrix of each substructure and connection coordinate transformation matrix and obtaining the FRF matrix of the substructure system
(3)
Calculating the connection receptance matrix of each connection substructure and obtaining the connection receptance matrix of the substructure system
(4)
Achieving the FRF matrix of the assembly structure by taking the obtained FRF matrix and connection receptance matrix of the substructure system into (15).


## 2.1. Synthesis Criterions of the Flexible Connection

Figure1 describes a substructure system divided from an assembly substructure with intricate connections that contains both flexible and rigid connections. Each connection of the substructure system is considered a substructure and is defined as a connection substructure (CS). A rigid connection essentially belongs to a special flexible connection in which stiffness tends towards infinity. Therefore, all the connections are considered flexible connections in the following derivation.Figure 1
Substructure system with intricate connections.We define two groups of connection coordinates in which a group corresponds to connection coordinatesb′, b¯, and b located on connection substructures, substructures, and assembly structure, respectively, and they are on a side of all the connections corresponding to another group, namely, c′, c¯, and c on another side. For a multidegree of freedom (DOF) connection substructure j, each DOF can be equivalent to a spring with equivalent stiffness-damping, as shown in Figure 2, in which each full line in the CS represents a DOF and the dotted line represents the remaining DOFs. The relation between response and excitation of the jth connection substructure is determined as follows:
(1)Fbj′Fcj′=Zj−Zj−ZjZjXbj′Xcj′,Zj=Kj+sCj,where j=1, 2,…,n, n is the sum of the connection substructures, Zj=zabo×o refers to the connection impedance matrix of the jth connection substructure, o denotes the DOF sum of the jth connection substructure, and zab denotes the connection impedance between ath and bth DOFs that denotes the force needed to act on the bth DOF when the ath DOF happens to be a unit displacement. When no couple exists between different connection DOFs, Zj corresponds to a diagonal matrix. Additionally, Kj and Cj represent the equivalent connection stiffness and damping matrixes, respectively, of the jth connection substructure; s=wi=2πfi, where i denotes the square root of −1 and w and f denote the vibration angular frequency and vibration frequency, respectively.Figure 2
Connection substructure between two substructures.Furthermore,Rj is defined as the connection receptance matrix of the jth connection substructure that is the inverse matrix of the connection impedance matrix and is calculated by (2). For a rigid connection, the connection receptance matrix is considered a zero matrix.
(2)Rj=Zj−1=Kj+sCj−1.Subsequently, the relation between response and excitation of thejth connection substructure is expressed by the connection receptance matrix as follows:
(3)Xbj′−Xcj′=RjFbj′=−RjFcj′.Similarly, the relation between response and excitation of all the connection substructures is expressed as follows:(4)Xb′−Xc′=RFb′=−RFc′,where R denotes the connection receptance matrix of the substructure system, and it is obtained by connection receptance matrixes of all the connection substructures in the diagonal form and expressed as follows:
(5)R=diagR1R2⋯Rn.


## 2.2. Formula Derivation

For the substructure system shown in Figure1, FRF matrixes of all the substructures are partitioned in internal coordinates and connection coordinates b¯ and c¯ that denote the relationship between response and excitation and are expressed as follows:
(6)Xa¯Xb¯Xc¯=H¯Fa¯Fb¯Fc¯=Ha¯a¯Ha¯b¯Ha¯c¯Hb¯a¯Hb¯b¯Hb¯c¯Hc¯a¯Hc¯b¯Hc¯c¯Fa¯Fb¯Fc¯.Accordingly, the relationship between response and excitation of the assembly structure is expressed as follows:(7)XaXbXc=HFaFbFc=HaaHabHacHbaHbbHbcHcaHcbHccFaFbFc,where subscripts a¯ and a represent the internal coordinates of the substructure system and the assembly structure, respectively.The basic conditions in the process of substructure synthesis are compatibility and equilibrium that are given as follows:(8)Xb¯Xc¯=Xb′Xc′=XbXc,(9)Fb¯Fc¯+Fb′Fc′=FbFc.In the substructure synthesis process, the response and excitation of the internal coordinates are not changed, and they are expressed as follows:(10)Xa¯=Xa,(11)Fa¯=Fa.By combining (4), (6), and (8), we obtain the following expression:
(12)RFb′=−RFc′=Hb¯a¯−Hc¯a¯Hb¯b¯−Hc¯b¯Hb¯c¯−Hc¯c¯TFa¯Fb¯Fc¯.The expression forFb′ and Fc′ is expressed as (13) by substituting (9) and (11) into (12) as follows:
(13)Fb′=−Fc′=R+Hb¯b¯+Hc¯c¯−Hc¯b¯−Hb¯c¯−1Hb¯a¯−Hc¯a¯Hb¯b¯−Hc¯b¯Hb¯c¯−Hc¯c¯TFaFbFc.Using (4), (8), (9), and (11), the relationship between response and excitation of the substructure system given by (6) is reexpressed as follows:
(14)XaXbXc=H¯FaFbFc+0−Fb′Fb′.By substituting (13) into (14) for Fb′, the FRF matrix of the assembly structure is obtained and expressed as (15) in which the FRF matrix of any assembly structure is calculated by the FRF and connection receptance matrixes of its substructures.
(15)H=H¯−Ha¯b¯−Ha¯c¯Hb¯b¯−Hb¯c¯Hb¯c¯−Hc¯c¯R+Hb¯b¯+Hc¯c¯−Hc¯b¯−Hb¯c¯−1Hb¯a¯−Hc¯a¯Hb¯b¯−Hc¯b¯Hb¯c¯−Hc¯c¯T.


## 2.3. Process of Substructure Synthesis

We consider a three-substructure system an example, and the process of substructure synthesis is demonstrated in this section. Figure3 shows the three-substructure system where substructure II connects to substructures I and III via connection substructures 1 and 2, respectively. Coordinates o1x1y1z1 (coordinate 1), o2x2y2z2 (coordinate 2), and o3x3y3z3 (coordinate 3) are the local coordinates of the three substructures, respectively, and are discordant from each other.Figure 3
Substructure system with three substructures.The FRF matrixes of the three substructures in their local coordinates are defined as follows:(16)HI=HaaIHab1IHb1aIHb1b1I,HII=HaaIIHac1IIHac2IIHc1aIIHc1c1IIHc1c2IIHc2aIIHc2c1IIHc2c2II,HIII=HaaIIIHab2IIIHb2aIIIHb2b2III.The FRF matrix of the substructure system is obtained by distributing FRF matrixes of the three substructures in internal coordinatea and connection coordinates b and c and expressed as follows:
(17)Due to the inconformity of the three coordinates, the connection coordinatesb and c are not in agreement. The FRF matrix in (17) is substituted into (15) only if the connection coordinates b and c are transformed into accordant coordinates. It is assumed that the connection coordinates of substructures I and III are accordant to those of substructure II after the transformation, and thus it is necessary to transform the connection coordinates of substructures I and III by a connection coordinate transformation matrix that is defined as follows:
(18)A=diagIuA1A2Iv,where I denotes the unit matrix, subscripts u and v refer to their order that is the sum of DOFs of all the internal and connection coordinates, and matrixes A1 and A2 represent the orientation cosine matrixes of all the connection coordinate DOFs of substructure II with those of substructures I and III, respectively.After the transformation via the connection coordinate transformation matrix, the FRF matrix of the substructure system is reexpressed as follows:(19)The connection receptance matrixes of connection substructures 1 and 2 are defined asR1 and R2, respectively, in coordinate 2. With respect to (5), the connection receptance matrix of the substructure system is expressed as follows:
(20)R=R100R2.By substituting (19) and (20) into (15), the FRF matrix of the assembly structure is achieved. The steps involved in the substructure synthesis process can be summarized as follows:
(1)
Dividing a complex assembly structure into several substructures depending on the analysis demand
(2)
Calculating the FRF matrix of each substructure and connection coordinate transformation matrix and obtaining the FRF matrix of the substructure system
(3)
Calculating the connection receptance matrix of each connection substructure and obtaining the connection receptance matrix of the substructure system
(4)
Achieving the FRF matrix of the assembly structure by taking the obtained FRF matrix and connection receptance matrix of the substructure system into (15).


## 3. Description and Simplification of the CMG Model

### 3.1. Description of the CMG Structure

The CMG structure mainly consists of the following four parts: flywheel, bearing, gimbal, and bracket [6], as shown in Figure 4. There are two groups of rotor-bearing systems in the CMG structure, namely, the flywheel-bearing system in which the flywheel rapidly spins to the gimbal and the gimbal-bearing system in which the gimbal slowly rotates to the base controlled by the gimbal servo system. The structural forms of both correspond to a rotor supported by a pair of face-to-face-mounted angular ball bearings. It is assumed centers of mass (COMs) of the flywheel, gimbal, and bracket are coincident at the initial state.Figure 4
Schematic diagram of the CMG structure.A few coordinates are defined for the purpose of convenience in modeling the CMG microvibrations. We defineo0x0y0z0 (coordinate 0) as an inertial coordinate fixed in the bracket in which the origin o0 is located at the COM of the bracket and the y0-axis is in line with the gimbal axis. We define ogxgygzg (coordinate g) as the gimbal-fixed coordinate that is a motion coordinate and rotates with the gimbal-bearing system. The origin og is located at the COM of the gimbal, the yg-axis is in line with the gimbal axis, and coordinate g initially coincides with coordinate 0.


### 3.2. Simplification of Rotor-Bearing Systems

For the two groups of rotor-bearing systems contained in the CMG structure, the flywheel, the gimbal, and the bracket are considered into six-DOF rigid bodies because the flywheel rotary speed is significantly lower than the critical point, and all the supporting bearings are equivalent to linear springs with radial and axial stiffness-damping [9]. Figure 5 shows the equivalent spring model of a rotor-bearing system in which kr and ka represent the radial and axial stiffness, respectively, of the equivalent linear spring, cr and ca represent the radial and axial damping, respectively, and d denotes the bearing supporting length that is the distance from COM of the rotor o to the bearing supporting point. All the radial and axial stiffness-damping are equivalent to the COM of the rotor, and the rotor-bearing system is simplified into a six-DOF mass-spring-damping vibration system in which stiffness and damping matrixes are defined as follows:
(21)Koxyz=diagKtrKto,Coxyz=diagCtrCto,where K, Ktr, and Kto represent the stiffness, translational stiffness, and torsional stiffness matrixes, respectively, of the simplified spring and C, Ctr, and Cto represent the damping, translational damping, and torsional damping matrixes, respectively.Figure 5
Equivalent spring model of the rotor-bearing system.Translational stiffness and damping matrixes of the vibration system are obtained by linear summation of the equivalent stiffness and damping of the two bearings, respectively, and given as follows:(22)Ktroxyz=diag2kr2kr2ka,Ctroxyz=diag2cr2cr2ca.Torsional stiffness matrix of the six-DOF vibration system is deduced by performing force analysis to the rotor-bearing system. The derivation of torsional stiffness around they-axis is considered an example to demonstrate this process. As depicted in Figure 6, it is assumed that the rotor-bearing system exhibits a small torsional deformation ε around the y-axis, and the two x-axis radial springs generate a translational transformation Δε in the opposite direction that is calculated as follows:
(23)Δε=dsinε.Figure 6
Equivalent model ofy-axis torsional springs.The reactive forces generated by the twox-axis radial springs are equal and opposite in direction. The value of the reactive forces is determined as follows:
(24)Fεoxyz=krΔε=krdsinε.The two reactive forces are equivalent to the origino as a torque around the y-axis that is obtained as follows:
(25)Tεoxyz=2Fεd=2krd2sinε.For the axis rotation of the flywheel and the gimbal about the two radial directions, the rotary angular are quite low, and thus the torsional deformationε tends towards zero. Subsequently, the equivalent torque is simplified as follows:
(26)Tεoxyz=2krd2ε.Thus, we obtain the equivalent torsional stiffness around they-axis as follows:
(27)ktoyoxyz=Tεε=2krd2.Similarly, the equivalent torsional stiffness around thex-axis is expressed as follows:
(28)ktoxoxyz=2krd2.The equivalent torsional stiffness around thez-axis is zero for the rotation axis of the rotor is the z-axis, and thus the equivalent torsional stiffness matrix is obtained as follows:
(29)Ktooxyz=diag2krd22krd20.Similarly, we obtain the equivalent torsional damping matrix of the system that is expressed as follows:(30)Ctooxyz=diag2crd22crd20.Finally, the equivalent stiffness and damping matrixes of the six-DOF vibration system are obtained as follows:(31)Koxyz=diag2kr2kr2ka2d2kr2d2kr0,Coxyz=diag2cr2cr2ca2d2cr2d2cr0.


### 3.3. Simplification of the CMG Model

Based on the simplification of the rotor-bearing system in the previous section, the two pairs of supporting bearings are equivalent to two linear springs with six DOFs, and the CMG model is simplified into an eighteen-DOF mass-spring-damping system as shown in Figure7.Figure 7
Simplified model of CMG.Based on (31), the equivalent stiffness and damping matrixes of the two equivalent linear springs are given as follows:
(32)K0ogxgygzg=diag2kr02kr02ka02d02kr02d02kr00,C0ogxgygzg=diag2cr02cr02ca02d02cr02d02cr00,(33)K1o0x0y0z0=diag2kr12ka12kr12d12kr102d12kr1,C1o0x0y0z0=diag2cr12ca12cr12d12cr102d12cr1,where subscripts 0 and 1 refer to the flywheel-bearing system and the gimbal-bearing system, respectively.The gimbal servo system is simplified into a torsional spring along the gimbal axis in which stiffness is defined as the servo dynamic stiffness [28]. Subsequently, the equivalent stiffness matrix in (33) is reexpressed as follows:
(34)K1o0x0y0z0=diag2kr12ka12kr12d12kr1ks2d12kr1,where ks represents the servo dynamic stiffness of the gimbal servo system in which the expression refers to [28].


## 3.1. Description of the CMG Structure

The CMG structure mainly consists of the following four parts: flywheel, bearing, gimbal, and bracket [6], as shown in Figure 4. There are two groups of rotor-bearing systems in the CMG structure, namely, the flywheel-bearing system in which the flywheel rapidly spins to the gimbal and the gimbal-bearing system in which the gimbal slowly rotates to the base controlled by the gimbal servo system. The structural forms of both correspond to a rotor supported by a pair of face-to-face-mounted angular ball bearings. It is assumed centers of mass (COMs) of the flywheel, gimbal, and bracket are coincident at the initial state.Figure 4
Schematic diagram of the CMG structure.A few coordinates are defined for the purpose of convenience in modeling the CMG microvibrations. We defineo0x0y0z0 (coordinate 0) as an inertial coordinate fixed in the bracket in which the origin o0 is located at the COM of the bracket and the y0-axis is in line with the gimbal axis. We define ogxgygzg (coordinate g) as the gimbal-fixed coordinate that is a motion coordinate and rotates with the gimbal-bearing system. The origin og is located at the COM of the gimbal, the yg-axis is in line with the gimbal axis, and coordinate g initially coincides with coordinate 0.


## 3.2. Simplification of Rotor-Bearing Systems

For the two groups of rotor-bearing systems contained in the CMG structure, the flywheel, the gimbal, and the bracket are considered into six-DOF rigid bodies because the flywheel rotary speed is significantly lower than the critical point, and all the supporting bearings are equivalent to linear springs with radial and axial stiffness-damping [9]. Figure 5 shows the equivalent spring model of a rotor-bearing system in which kr and ka represent the radial and axial stiffness, respectively, of the equivalent linear spring, cr and ca represent the radial and axial damping, respectively, and d denotes the bearing supporting length that is the distance from COM of the rotor o to the bearing supporting point. All the radial and axial stiffness-damping are equivalent to the COM of the rotor, and the rotor-bearing system is simplified into a six-DOF mass-spring-damping vibration system in which stiffness and damping matrixes are defined as follows:
(21)Koxyz=diagKtrKto,Coxyz=diagCtrCto,where K, Ktr, and Kto represent the stiffness, translational stiffness, and torsional stiffness matrixes, respectively, of the simplified spring and C, Ctr, and Cto represent the damping, translational damping, and torsional damping matrixes, respectively.Figure 5
Equivalent spring model of the rotor-bearing system.Translational stiffness and damping matrixes of the vibration system are obtained by linear summation of the equivalent stiffness and damping of the two bearings, respectively, and given as follows:(22)Ktroxyz=diag2kr2kr2ka,Ctroxyz=diag2cr2cr2ca.Torsional stiffness matrix of the six-DOF vibration system is deduced by performing force analysis to the rotor-bearing system. The derivation of torsional stiffness around they-axis is considered an example to demonstrate this process. As depicted in Figure 6, it is assumed that the rotor-bearing system exhibits a small torsional deformation ε around the y-axis, and the two x-axis radial springs generate a translational transformation Δε in the opposite direction that is calculated as follows:
(23)Δε=dsinε.Figure 6
Equivalent model ofy-axis torsional springs.The reactive forces generated by the twox-axis radial springs are equal and opposite in direction. The value of the reactive forces is determined as follows:
(24)Fεoxyz=krΔε=krdsinε.The two reactive forces are equivalent to the origino as a torque around the y-axis that is obtained as follows:
(25)Tεoxyz=2Fεd=2krd2sinε.For the axis rotation of the flywheel and the gimbal about the two radial directions, the rotary angular are quite low, and thus the torsional deformationε tends towards zero. Subsequently, the equivalent torque is simplified as follows:
(26)Tεoxyz=2krd2ε.Thus, we obtain the equivalent torsional stiffness around they-axis as follows:
(27)ktoyoxyz=Tεε=2krd2.Similarly, the equivalent torsional stiffness around thex-axis is expressed as follows:
(28)ktoxoxyz=2krd2.The equivalent torsional stiffness around thez-axis is zero for the rotation axis of the rotor is the z-axis, and thus the equivalent torsional stiffness matrix is obtained as follows:
(29)Ktooxyz=diag2krd22krd20.Similarly, we obtain the equivalent torsional damping matrix of the system that is expressed as follows:(30)Ctooxyz=diag2crd22crd20.Finally, the equivalent stiffness and damping matrixes of the six-DOF vibration system are obtained as follows:(31)Koxyz=diag2kr2kr2ka2d2kr2d2kr0,Coxyz=diag2cr2cr2ca2d2cr2d2cr0.


## 3.3. Simplification of the CMG Model

Based on the simplification of the rotor-bearing system in the previous section, the two pairs of supporting bearings are equivalent to two linear springs with six DOFs, and the CMG model is simplified into an eighteen-DOF mass-spring-damping system as shown in Figure7.Figure 7
Simplified model of CMG.Based on (31), the equivalent stiffness and damping matrixes of the two equivalent linear springs are given as follows:
(32)K0ogxgygzg=diag2kr02kr02ka02d02kr02d02kr00,C0ogxgygzg=diag2cr02cr02ca02d02cr02d02cr00,(33)K1o0x0y0z0=diag2kr12ka12kr12d12kr102d12kr1,C1o0x0y0z0=diag2cr12ca12cr12d12cr102d12cr1,where subscripts 0 and 1 refer to the flywheel-bearing system and the gimbal-bearing system, respectively.The gimbal servo system is simplified into a torsional spring along the gimbal axis in which stiffness is defined as the servo dynamic stiffness [28]. Subsequently, the equivalent stiffness matrix in (33) is reexpressed as follows:
(34)K1o0x0y0z0=diag2kr12ka12kr12d12kr1ks2d12kr1,where ks represents the servo dynamic stiffness of the gimbal servo system in which the expression refers to [28].


## 4. Dynamic Responses of CMG on the Flexible Interface

### 4.1. Division of the Substructures

We consider an aluminum honeycomb sandwich palate (AHSP) a flexible interface, and the CMG-flexible interface coupling system is shown in Figure8 in which the CMG is hard-mounted on the AHSP. Given the complicated kinematics of the CMG microvibration system, the CMG structure is divided into the following two parts: the flywheel-gimbal part and bracket part in which it is more convenient to independently develop equations of motion (EOMs), and this is followed by the complete CMG. The flexible interface is separated into a substructure. Therefore, the coupling system is divided into the following three substructures: the flywheel-gimbal substructure (substructure I), bracket substructure (substructure II), and flexible interface substructure (substructure III) as shown in Figure 9, and their local coordinates are coordinate g, coordinate 0, and coordinate oaxayaza (coordinate a), respectively. Coordinate a is defined as the interface between the CMG and the flexible mounting surface, and the origin oa is located at the center of the installation surface. The connection bearings between the gimbal and the bracket are considered a six-DOF flexible connection (connection substructure 1), and the rigid connection between the bracket and the AHSP is considered a rigid connection (connection substructure 2).Figure 8
CMG-flexible interface coupling system.Figure 9
Substructure division of the coupling system.


### 4.2. Substructure I

Substructure I consists of the flywheel and the gimbal connected by a pair of angular ball bearings, and this is simplified into a twelve-DOF mass-spring-damping system based on the previous simplification. This part develops the EOM of substructure I and uses it to calculate the FRF matrix.


#### 4.2.1. Modeling Method

The Lagrange energy method is used to develop the dynamic model of substructure I. The general form of the Lagrange equation is given as follows:(35)ddt∂Tq,q̇∂q̇−∂Tq,q̇∂q+∂Vq,q̇∂q+∂Dq,q̇∂q̇=p,where q is the vector of the generalized freedom, the terms T, V, and D denote the kinetic energy, potentials, and dissipations of the system, respectively, and p represents the generalized force vector.


#### 4.2.2. Kinetic Energy

The displacement vectors of the COM of the flywheel and the gimbal in coordinateg are defined as follows:
(36)qfogxgygzg=qftTqfrTT=xfyfzfαfβfγfT,qgogxgygzg=qgtTqgrTT=xgygzgαgβgγgT,where qft and qgt are the translational displacement vectors of the COM of the flywheel and the gimbal, respectively, and qfr and qgr are the angular displacement vectors, respectively.The velocity vectors are derived from differential of the displacement vectors and expressed as follows:(37)q̇fogxgygzg=q̇ftTq̇frTT=ẋfẏfżfα̇fβ̇fγ̇fT,q̇gogxgygzg=q̇gtTq̇grTT=ẋgẏgżgα̇gβ̇gγ̇gT,where q̇ft and q̇gt are the translational velocity vectors of the COM of the flywheel and the gimbal, respectively, in coordinate g and q̇fr and q̇gr are the angular velocity vectors, respectively.As shown in Figure10, generalized angular motions are used to define the rigid body rotations of the flywheel and the gimbal that are shown by the Euler angles. The first rotation, β, is about the yg-axis of the inertial coordinate, namely, coordinate g, and it defines the intermediate reference coordinate, namely, oβxβyβzβ (coordinate β). The next rotation, α, that is about the xβ-axis defines the rocking frame, namely, oαxαyαzα (coordinate α). The final rotation, γ, is about the zα-axis and defines the final body-fixed coordinate, namely, oγxγyγzγ (coordinate γ).Figure 10
Euler angles of the rocking motion of the flywheel and the gimbal.After the generalized angular motion, the angular velocity vectors of the COM of the flywheel and the gimbal in the rocking frame,oαxαyαzα, are expressed as follows:
(38)q̇froαxαyαzα=α̇fβ̇fcosβfγ̇f−β̇fsinβfT,q̇groαxαyαzα=α̇gβ̇gcosβgγ̇g−β̇gsinβgT.It is assumed that the flywheel spinning speed is a constantΩ, and the angular velocity vector q̇fr is reexpressed as follows:
(39)q̇froαxαyαzα=α̇fβ̇fcosβfΩ−β̇fsinβfT.The mass matrixes of the flywheel and the gimbal are given as follows:(40)Mf=diagmfmfmfJfrJfrJfz,Mg=diagmgmgmgJgrJgrJgz,where mf and mg are the mass of the flywheel and the gimbal, respectively, Jfr and Jgr denote the radial inertia, and Jfz and Jgz denote the polar inertia.The total kinetic energyT of the system is determined as follows:
(41)T=12q̇fTMfq̇f+12q̇gTMgq̇g.Finally, the total kinetic energy of the system is obtained by substituting (36), (37), (38), (39), and (40) into (41) and expressed as follows:
(42)T=12Mgẋg2+ẏg2+żg2+12Jgrα̇g2+β̇gcosβg2+12Jgzγ̇g−β̇gsinβg2+12Mfẋf2+ẏf2+żf2+12Jfrα̇f2+β̇fcosβf2+12JfzΩ−β̇fsinβf2.


#### 4.2.3. Potentials and Dissipations

The stiffness and damping matrixes of the vibration system are given by (32). The elastic displacement vector of the system is obtained by subtracting the displacement vector of the COM of the gimbal from that of the flywheel as follows:
(43)qgfogxgygzg=qg−qf.The potentials of the system are calculated as follows:(44)V=12qgfTK1qgf.Thus, we obtain the potentials of the system by combining (32), (43), and (44) that are expressed as follows:
(45)V=kr1xf−xg2+yf−yg2+ka1zf−zg2+d12kr1αf−αg2+βf−βg2.The dissipation velocity vector of the system is calculated by subtracting the velocity vector of the COM of the gimbal from that of the flywheel as follows:(46)q̇gfogxgygzg=q̇g−q̇f.The dissipations of the system are determined as follows:(47)D=12q̇gfTC1q̇gf.Hence, (32), (46), and (47) are used to obtain the dissipations of the system as follows:
(48)D=cr1ẋf−ẋg2+ẏf−ẏg2+ca1żf−żg2+d12cr1α̇f−α̇g2+β̇f−β̇g2.


#### 4.2.4. Generalized Force

The generalized external force vector of the system is defined as follows:(49)pIogxgygzg=ffTfgTT,where ff and fg represent the generalized external force vectors acting on the COM of the flywheel and the gimbal, respectively.The generalized external force vectors acting on the flywheel are mainly attributed to the imbalance features caused by the flywheel mass imbalance. There are two types of imbalance features: static and dynamic imbalances, and these are shown in Figure11.Figure 11
The static and dynamic imbalances of the flywheel.Static imbalance is caused by the offset of the flywheel COM from the rotation axis. The dynamic imbalance results from the angular misalignment of the principle axis of the flywheel and the rotation axis. When the wheel rotates, the static and dynamic imbalances cause imbalance forces and torques to the flywheel that are equivalent to a force vector and a torque vector, respectively, acting on the COM of the flywheel. The equivalent force and torque vectors are expressed as follows:(50)fsogxgygzg=UsΩ2sinΩt+φsUsΩ2cosΩt+φs0T,fdogxgygzg=UdΩ2sinΩt+φdUdΩ2cosΩt+φd0T,where Us=msrs and Ud=mdrdhd denote the static and dynamic mass imbalances and φs and φd denote the initial phases, respectively.Subsequently, the generalized external force vector acting on the COM of the flywheel is expressed as follows:(51)ffogxgygzg=fsTfdTT.


#### 4.2.5. FRF Matrix of Substructure I

By substituting (42), (45), (48), and (49) into (35), we obtain the EOM of substructure I in the time domain. Terms with order exceeding O(qi2) are ignored, and the EOM is expressed in a matrix form as follows:
(52)MIq¨I+CI+GIq̇I+KIqI=pI,where MI, CI, and KI represent the generalized mass, damping, and stiffness matrixes, respectively, of the system and GI denotes the centrifugal term matrix, and all the expressions are expressed as follows:
(53)MI=Mf00Mg,KI=K0−K0−K0K0,CI=C0−C0−C0C0,GI=Gf000,where Gf is the centrifugal matrix of the flywheel that is calculated as follows:
(54)Gf=000G1,G1=0JfzΩ0−JfzΩ00000.The generalized coordinate vector of the vibration system is expressed as follows:(55)qIogxgygzg=qfTqgTT=xfyfzfαfβfγfxgygzgαgβgγgT,where qf and qg are the displacement vectors of the COM of the flywheel and the gimbal that denote the internal and connection coordinate vectors of substructure I, respectively.The EOM in (52) is transformed into the frequency domain form via the Fourier transform (FT) that is expressed as follows:
(56)s2MIQIs+sCI+GIQIs+KIQIs=PIs.We arrange (56) as follows:
(57)ZIQIs=PIs,ZI=s2MI+sCI+GI+KI,where ZI is the impedance matrix of the vibration system and is generally invertible, and the FRF matrix of substructure I is derived from the impedance matrix as follows:
(58)HI=ZI−1=HffHfgHgfHgg.


### 4.3. Substructure II

Substructure II is the bracket substructure that is considered a six-DOF rigid body. Two connection coordinates are contained in substructure II in which a coordinate is located on the COM of the bracket connecting to connection substructure 1 and is defined asc1 and another coordinate is located on the center of the bracket undersurface connecting to connection substructure 2 and is defined as c2.Substructure II is a six-DOF rigid body, and thus the response and excitation vectors between the two connection coordinates are equivalent to each other. The equivalent relations are determined by the following formulae:(59)Fc1II=B1Fc2II,Xc2II=B1TXc1II,B1=I0B2I,B2=0h0h00000,where h is the distance between the two connection coordinates c1 and c2 and matrix B1 and its transpose B1T are transform matrixes of the response and excitation vectors between the two connection coordinates, respectively.The response and excitation of coordinatec1 follow this rule and are expressed as follows:
(60)Xc1II=HcFc1II,Hc=s2Mb−1,Mb=diagmbmbmbJbrJbrJbz,where mb is the mass of the bracket and Jbr and Jbz are the radial inertia and polar inertia, respectively.The transitive relation between response and excitation for the two connection coordinates is derived from (59) and (60) that is given as follows:
(61)Xc1IIXc2II=HcHcB1B1THcB1THcB1Fc1IIFc2II.Therefore, we obtain the FRF matrix of substructureII as follows:
(62)HII=HcHcB1B1THcB1THcB1.


### 4.4. Substructure III

Substructure III is the flexible interface that is demonstrated by an AHSP. Substructure III only considers a connection coordinate (coordinateb2) for the study and mainly investigates the dynamic responses of the CMG mounting point. The FRF matrix of substructure III is obtained by its finite element model (FEM).The structure of the AHSP consists of two layers of aluminum skins on the two sides with honeycomb core in the middle. The honeycomb core is equivalent to homogeneous orthotropic material, and thus the AHSP is considered a laminate with three-layer plates. The size of the AHSP is 1 m × 0.5 m × 0.03 m, and the thickness of a single-layer aluminum skin is 0.3 mm. Additionally,E, G, ρ, and υ denote Young’s modulus, shear modulus, density, and Poisson ratio, respectively. The material properties of the aluminum skins and the honeycomb core equivalent model are shown in Tables 1 and 2, respectively, and are derived from reference [29].Table 1
Material properties of the aluminum alloy.
Young’s modulus (E)
Density (ρ)
Poisson ratio (υ)
71 GPa
2700 kg/m3
0.33Table 2
Equivalent material properties of the honeycomb core.
Symbol
Values
Unit
Exx
0.0354
MPa
Eyy
0.0354
MPa
Ezz
655.87
MPa
Gxy
0.0266
MPa
Gyz
141.12
MPa
Gzx
92.463
MPa
ρ
24.94
kg/m3
υ12
0.999856
—
υ23
0
—
υ31
0
—Furthermore, MSC PATRAN is used to establish the FEM of the AHSP, as shown in Figure12. In the FEM, a laminate element is applied to mesh the AHSP. The equivalent material properties of the honeycomb core are simulated using a 2D orthotropic material. The nodes on the four sides of the FEM are imposed by six-DOF displacement constraints, and the separation distance between two nodes is 0.1 m. The FRF matrix of substructure III Hb2 is obtained via harmonic response analysis of the FEM by acting six-DOF unit force on the CMG mounting point independently.Figure 12
Schematic diagram of FEM of the AHSP.


### 4.5. FRF Matrix of the Coupling System

With respect to the substructure system, connection coordinates on the two sides of connection substructure 1 are not in agreement. In the process of providing control moments, the gimbal rotates corresponding to the bracket around the shaft, and thus coordinateg rotates to coordinate 0 around the y-axis as shown in Figure 13 in which θ is the rotation angle. Table 3 lists the angle relation between coordinates g and 0 in which the orientation cosine matrix is obtained as follows:
(63)A1=A300A3,A3=cosθ0cos0.5π−θ010cos0.5π+θ0cosθ.Figure 13
Angle relation between coordinatesg and 0.Table 3
Angle relation between the coordinatesg and 0.
Axis
xg
yg
zg
x0
θ
0.5π
0.5π − θ
y0
0.5π
0
0.5π
z0
0.5π + θ
0.5π
θConnection coordinates on the two sides of connection substructure 2 are accordant, and thus the orientation cosine matrix between coordinatesa and 0 is the unit matrix. The sum of the DOFs of all the internal and connection coordinates of the substructure system is six and twelve, respectively. Subsequently, the connection coordinate transformation matrix of the substructure system is obtained and expressed as follows:
(64)AF=diagI6A1I6I12.Finally, the FRF matrix of the substructure system is achieved by substituting (58), (62), and (64) into (19) and is expressed as follows:
(65)Moreover, (2), (33), and (34) are used, and the connection receptance matrix of connection substructure 1 is expressed as follows:
(66)R1o0x0y0z0=Z1−1=K1+sC1−1.Connection receptance of connection substructure 2 tends towards zero because it is a rigid connection. Hence, the connection receptance matrix of the substructure system is expressed as follows:(67)RF=Z1−1000.Finally, the FRF matrix of the coupling systemHF is achieved by substituting (65) and (67) into (15).


### 4.6. Dynamic Responses of the Coupling System

For the coupling system, only external excitation acts on the flywheel as discussed in Section4.2.4. and shown in (50). Thus, the external excitation vector of the coupling system is obtained as follows:
(68)fF=fsTfdT0TT.We apply the FT to (68), and the external excitation vector is transformed into the frequency domain form that is expressed as follows:
(69)FFs=FsTsFdTs0TT.Displacement response of the coupling system is obtained by using its FRF matrix to multiply the external excitation vector as follows:(70)XF=HFFF.The displacement response is multiplied bys2, and we obtain the acceleration response of the coupling system as follows:
(71)X¨F=HFFFs2.Reaction forces and moments between CMG and the flexible interface are achieved by the impedance matrix of the CMG mounting point multiplied by its displacement response vector as follows:(72)Fm=Hb1−1Xm,where Xm (a part of XF) is the displacement response vector of the CMG mounting point.


## 4.1. Division of the Substructures

We consider an aluminum honeycomb sandwich palate (AHSP) a flexible interface, and the CMG-flexible interface coupling system is shown in Figure8 in which the CMG is hard-mounted on the AHSP. Given the complicated kinematics of the CMG microvibration system, the CMG structure is divided into the following two parts: the flywheel-gimbal part and bracket part in which it is more convenient to independently develop equations of motion (EOMs), and this is followed by the complete CMG. The flexible interface is separated into a substructure. Therefore, the coupling system is divided into the following three substructures: the flywheel-gimbal substructure (substructure I), bracket substructure (substructure II), and flexible interface substructure (substructure III) as shown in Figure 9, and their local coordinates are coordinate g, coordinate 0, and coordinate oaxayaza (coordinate a), respectively. Coordinate a is defined as the interface between the CMG and the flexible mounting surface, and the origin oa is located at the center of the installation surface. The connection bearings between the gimbal and the bracket are considered a six-DOF flexible connection (connection substructure 1), and the rigid connection between the bracket and the AHSP is considered a rigid connection (connection substructure 2).Figure 8
CMG-flexible interface coupling system.Figure 9
Substructure division of the coupling system.


## 4.2. Substructure I

Substructure I consists of the flywheel and the gimbal connected by a pair of angular ball bearings, and this is simplified into a twelve-DOF mass-spring-damping system based on the previous simplification. This part develops the EOM of substructure I and uses it to calculate the FRF matrix.


### 4.2.1. Modeling Method

The Lagrange energy method is used to develop the dynamic model of substructure I. The general form of the Lagrange equation is given as follows:(35)ddt∂Tq,q̇∂q̇−∂Tq,q̇∂q+∂Vq,q̇∂q+∂Dq,q̇∂q̇=p,where q is the vector of the generalized freedom, the terms T, V, and D denote the kinetic energy, potentials, and dissipations of the system, respectively, and p represents the generalized force vector.


### 4.2.2. Kinetic Energy

The displacement vectors of the COM of the flywheel and the gimbal in coordinateg are defined as follows:
(36)qfogxgygzg=qftTqfrTT=xfyfzfαfβfγfT,qgogxgygzg=qgtTqgrTT=xgygzgαgβgγgT,where qft and qgt are the translational displacement vectors of the COM of the flywheel and the gimbal, respectively, and qfr and qgr are the angular displacement vectors, respectively.The velocity vectors are derived from differential of the displacement vectors and expressed as follows:(37)q̇fogxgygzg=q̇ftTq̇frTT=ẋfẏfżfα̇fβ̇fγ̇fT,q̇gogxgygzg=q̇gtTq̇grTT=ẋgẏgżgα̇gβ̇gγ̇gT,where q̇ft and q̇gt are the translational velocity vectors of the COM of the flywheel and the gimbal, respectively, in coordinate g and q̇fr and q̇gr are the angular velocity vectors, respectively.As shown in Figure10, generalized angular motions are used to define the rigid body rotations of the flywheel and the gimbal that are shown by the Euler angles. The first rotation, β, is about the yg-axis of the inertial coordinate, namely, coordinate g, and it defines the intermediate reference coordinate, namely, oβxβyβzβ (coordinate β). The next rotation, α, that is about the xβ-axis defines the rocking frame, namely, oαxαyαzα (coordinate α). The final rotation, γ, is about the zα-axis and defines the final body-fixed coordinate, namely, oγxγyγzγ (coordinate γ).Figure 10
Euler angles of the rocking motion of the flywheel and the gimbal.After the generalized angular motion, the angular velocity vectors of the COM of the flywheel and the gimbal in the rocking frame,oαxαyαzα, are expressed as follows:
(38)q̇froαxαyαzα=α̇fβ̇fcosβfγ̇f−β̇fsinβfT,q̇groαxαyαzα=α̇gβ̇gcosβgγ̇g−β̇gsinβgT.It is assumed that the flywheel spinning speed is a constantΩ, and the angular velocity vector q̇fr is reexpressed as follows:
(39)q̇froαxαyαzα=α̇fβ̇fcosβfΩ−β̇fsinβfT.The mass matrixes of the flywheel and the gimbal are given as follows:(40)Mf=diagmfmfmfJfrJfrJfz,Mg=diagmgmgmgJgrJgrJgz,where mf and mg are the mass of the flywheel and the gimbal, respectively, Jfr and Jgr denote the radial inertia, and Jfz and Jgz denote the polar inertia.The total kinetic energyT of the system is determined as follows:
(41)T=12q̇fTMfq̇f+12q̇gTMgq̇g.Finally, the total kinetic energy of the system is obtained by substituting (36), (37), (38), (39), and (40) into (41) and expressed as follows:
(42)T=12Mgẋg2+ẏg2+żg2+12Jgrα̇g2+β̇gcosβg2+12Jgzγ̇g−β̇gsinβg2+12Mfẋf2+ẏf2+żf2+12Jfrα̇f2+β̇fcosβf2+12JfzΩ−β̇fsinβf2.


### 4.2.3. Potentials and Dissipations

The stiffness and damping matrixes of the vibration system are given by (32). The elastic displacement vector of the system is obtained by subtracting the displacement vector of the COM of the gimbal from that of the flywheel as follows:
(43)qgfogxgygzg=qg−qf.The potentials of the system are calculated as follows:(44)V=12qgfTK1qgf.Thus, we obtain the potentials of the system by combining (32), (43), and (44) that are expressed as follows:
(45)V=kr1xf−xg2+yf−yg2+ka1zf−zg2+d12kr1αf−αg2+βf−βg2.The dissipation velocity vector of the system is calculated by subtracting the velocity vector of the COM of the gimbal from that of the flywheel as follows:(46)q̇gfogxgygzg=q̇g−q̇f.The dissipations of the system are determined as follows:(47)D=12q̇gfTC1q̇gf.Hence, (32), (46), and (47) are used to obtain the dissipations of the system as follows:
(48)D=cr1ẋf−ẋg2+ẏf−ẏg2+ca1żf−żg2+d12cr1α̇f−α̇g2+β̇f−β̇g2.


### 4.2.4. Generalized Force

The generalized external force vector of the system is defined as follows:(49)pIogxgygzg=ffTfgTT,where ff and fg represent the generalized external force vectors acting on the COM of the flywheel and the gimbal, respectively.The generalized external force vectors acting on the flywheel are mainly attributed to the imbalance features caused by the flywheel mass imbalance. There are two types of imbalance features: static and dynamic imbalances, and these are shown in Figure11.Figure 11
The static and dynamic imbalances of the flywheel.Static imbalance is caused by the offset of the flywheel COM from the rotation axis. The dynamic imbalance results from the angular misalignment of the principle axis of the flywheel and the rotation axis. When the wheel rotates, the static and dynamic imbalances cause imbalance forces and torques to the flywheel that are equivalent to a force vector and a torque vector, respectively, acting on the COM of the flywheel. The equivalent force and torque vectors are expressed as follows:(50)fsogxgygzg=UsΩ2sinΩt+φsUsΩ2cosΩt+φs0T,fdogxgygzg=UdΩ2sinΩt+φdUdΩ2cosΩt+φd0T,where Us=msrs and Ud=mdrdhd denote the static and dynamic mass imbalances and φs and φd denote the initial phases, respectively.Subsequently, the generalized external force vector acting on the COM of the flywheel is expressed as follows:(51)ffogxgygzg=fsTfdTT.


### 4.2.5. FRF Matrix of Substructure I

By substituting (42), (45), (48), and (49) into (35), we obtain the EOM of substructure I in the time domain. Terms with order exceeding O(qi2) are ignored, and the EOM is expressed in a matrix form as follows:
(52)MIq¨I+CI+GIq̇I+KIqI=pI,where MI, CI, and KI represent the generalized mass, damping, and stiffness matrixes, respectively, of the system and GI denotes the centrifugal term matrix, and all the expressions are expressed as follows:
(53)MI=Mf00Mg,KI=K0−K0−K0K0,CI=C0−C0−C0C0,GI=Gf000,where Gf is the centrifugal matrix of the flywheel that is calculated as follows:
(54)Gf=000G1,G1=0JfzΩ0−JfzΩ00000.The generalized coordinate vector of the vibration system is expressed as follows:(55)qIogxgygzg=qfTqgTT=xfyfzfαfβfγfxgygzgαgβgγgT,where qf and qg are the displacement vectors of the COM of the flywheel and the gimbal that denote the internal and connection coordinate vectors of substructure I, respectively.The EOM in (52) is transformed into the frequency domain form via the Fourier transform (FT) that is expressed as follows:
(56)s2MIQIs+sCI+GIQIs+KIQIs=PIs.We arrange (56) as follows:
(57)ZIQIs=PIs,ZI=s2MI+sCI+GI+KI,where ZI is the impedance matrix of the vibration system and is generally invertible, and the FRF matrix of substructure I is derived from the impedance matrix as follows:
(58)HI=ZI−1=HffHfgHgfHgg.


## 4.2.1. Modeling Method

The Lagrange energy method is used to develop the dynamic model of substructure I. The general form of the Lagrange equation is given as follows:(35)ddt∂Tq,q̇∂q̇−∂Tq,q̇∂q+∂Vq,q̇∂q+∂Dq,q̇∂q̇=p,where q is the vector of the generalized freedom, the terms T, V, and D denote the kinetic energy, potentials, and dissipations of the system, respectively, and p represents the generalized force vector.


## 4.2.2. Kinetic Energy

The displacement vectors of the COM of the flywheel and the gimbal in coordinateg are defined as follows:
(36)qfogxgygzg=qftTqfrTT=xfyfzfαfβfγfT,qgogxgygzg=qgtTqgrTT=xgygzgαgβgγgT,where qft and qgt are the translational displacement vectors of the COM of the flywheel and the gimbal, respectively, and qfr and qgr are the angular displacement vectors, respectively.The velocity vectors are derived from differential of the displacement vectors and expressed as follows:(37)q̇fogxgygzg=q̇ftTq̇frTT=ẋfẏfżfα̇fβ̇fγ̇fT,q̇gogxgygzg=q̇gtTq̇grTT=ẋgẏgżgα̇gβ̇gγ̇gT,where q̇ft and q̇gt are the translational velocity vectors of the COM of the flywheel and the gimbal, respectively, in coordinate g and q̇fr and q̇gr are the angular velocity vectors, respectively.As shown in Figure10, generalized angular motions are used to define the rigid body rotations of the flywheel and the gimbal that are shown by the Euler angles. The first rotation, β, is about the yg-axis of the inertial coordinate, namely, coordinate g, and it defines the intermediate reference coordinate, namely, oβxβyβzβ (coordinate β). The next rotation, α, that is about the xβ-axis defines the rocking frame, namely, oαxαyαzα (coordinate α). The final rotation, γ, is about the zα-axis and defines the final body-fixed coordinate, namely, oγxγyγzγ (coordinate γ).Figure 10
Euler angles of the rocking motion of the flywheel and the gimbal.After the generalized angular motion, the angular velocity vectors of the COM of the flywheel and the gimbal in the rocking frame,oαxαyαzα, are expressed as follows:
(38)q̇froαxαyαzα=α̇fβ̇fcosβfγ̇f−β̇fsinβfT,q̇groαxαyαzα=α̇gβ̇gcosβgγ̇g−β̇gsinβgT.It is assumed that the flywheel spinning speed is a constantΩ, and the angular velocity vector q̇fr is reexpressed as follows:
(39)q̇froαxαyαzα=α̇fβ̇fcosβfΩ−β̇fsinβfT.The mass matrixes of the flywheel and the gimbal are given as follows:(40)Mf=diagmfmfmfJfrJfrJfz,Mg=diagmgmgmgJgrJgrJgz,where mf and mg are the mass of the flywheel and the gimbal, respectively, Jfr and Jgr denote the radial inertia, and Jfz and Jgz denote the polar inertia.The total kinetic energyT of the system is determined as follows:
(41)T=12q̇fTMfq̇f+12q̇gTMgq̇g.Finally, the total kinetic energy of the system is obtained by substituting (36), (37), (38), (39), and (40) into (41) and expressed as follows:
(42)T=12Mgẋg2+ẏg2+żg2+12Jgrα̇g2+β̇gcosβg2+12Jgzγ̇g−β̇gsinβg2+12Mfẋf2+ẏf2+żf2+12Jfrα̇f2+β̇fcosβf2+12JfzΩ−β̇fsinβf2.


## 4.2.3. Potentials and Dissipations

The stiffness and damping matrixes of the vibration system are given by (32). The elastic displacement vector of the system is obtained by subtracting the displacement vector of the COM of the gimbal from that of the flywheel as follows:
(43)qgfogxgygzg=qg−qf.The potentials of the system are calculated as follows:(44)V=12qgfTK1qgf.Thus, we obtain the potentials of the system by combining (32), (43), and (44) that are expressed as follows:
(45)V=kr1xf−xg2+yf−yg2+ka1zf−zg2+d12kr1αf−αg2+βf−βg2.The dissipation velocity vector of the system is calculated by subtracting the velocity vector of the COM of the gimbal from that of the flywheel as follows:(46)q̇gfogxgygzg=q̇g−q̇f.The dissipations of the system are determined as follows:(47)D=12q̇gfTC1q̇gf.Hence, (32), (46), and (47) are used to obtain the dissipations of the system as follows:
(48)D=cr1ẋf−ẋg2+ẏf−ẏg2+ca1żf−żg2+d12cr1α̇f−α̇g2+β̇f−β̇g2.


## 4.2.4. Generalized Force

The generalized external force vector of the system is defined as follows:(49)pIogxgygzg=ffTfgTT,where ff and fg represent the generalized external force vectors acting on the COM of the flywheel and the gimbal, respectively.The generalized external force vectors acting on the flywheel are mainly attributed to the imbalance features caused by the flywheel mass imbalance. There are two types of imbalance features: static and dynamic imbalances, and these are shown in Figure11.Figure 11
The static and dynamic imbalances of the flywheel.Static imbalance is caused by the offset of the flywheel COM from the rotation axis. The dynamic imbalance results from the angular misalignment of the principle axis of the flywheel and the rotation axis. When the wheel rotates, the static and dynamic imbalances cause imbalance forces and torques to the flywheel that are equivalent to a force vector and a torque vector, respectively, acting on the COM of the flywheel. The equivalent force and torque vectors are expressed as follows:(50)fsogxgygzg=UsΩ2sinΩt+φsUsΩ2cosΩt+φs0T,fdogxgygzg=UdΩ2sinΩt+φdUdΩ2cosΩt+φd0T,where Us=msrs and Ud=mdrdhd denote the static and dynamic mass imbalances and φs and φd denote the initial phases, respectively.Subsequently, the generalized external force vector acting on the COM of the flywheel is expressed as follows:(51)ffogxgygzg=fsTfdTT.


## 4.2.5. FRF Matrix of Substructure I

By substituting (42), (45), (48), and (49) into (35), we obtain the EOM of substructure I in the time domain. Terms with order exceeding O(qi2) are ignored, and the EOM is expressed in a matrix form as follows:
(52)MIq¨I+CI+GIq̇I+KIqI=pI,where MI, CI, and KI represent the generalized mass, damping, and stiffness matrixes, respectively, of the system and GI denotes the centrifugal term matrix, and all the expressions are expressed as follows:
(53)MI=Mf00Mg,KI=K0−K0−K0K0,CI=C0−C0−C0C0,GI=Gf000,where Gf is the centrifugal matrix of the flywheel that is calculated as follows:
(54)Gf=000G1,G1=0JfzΩ0−JfzΩ00000.The generalized coordinate vector of the vibration system is expressed as follows:(55)qIogxgygzg=qfTqgTT=xfyfzfαfβfγfxgygzgαgβgγgT,where qf and qg are the displacement vectors of the COM of the flywheel and the gimbal that denote the internal and connection coordinate vectors of substructure I, respectively.The EOM in (52) is transformed into the frequency domain form via the Fourier transform (FT) that is expressed as follows:
(56)s2MIQIs+sCI+GIQIs+KIQIs=PIs.We arrange (56) as follows:
(57)ZIQIs=PIs,ZI=s2MI+sCI+GI+KI,where ZI is the impedance matrix of the vibration system and is generally invertible, and the FRF matrix of substructure I is derived from the impedance matrix as follows:
(58)HI=ZI−1=HffHfgHgfHgg.


## 4.3. Substructure II

Substructure II is the bracket substructure that is considered a six-DOF rigid body. Two connection coordinates are contained in substructure II in which a coordinate is located on the COM of the bracket connecting to connection substructure 1 and is defined asc1 and another coordinate is located on the center of the bracket undersurface connecting to connection substructure 2 and is defined as c2.Substructure II is a six-DOF rigid body, and thus the response and excitation vectors between the two connection coordinates are equivalent to each other. The equivalent relations are determined by the following formulae:(59)Fc1II=B1Fc2II,Xc2II=B1TXc1II,B1=I0B2I,B2=0h0h00000,where h is the distance between the two connection coordinates c1 and c2 and matrix B1 and its transpose B1T are transform matrixes of the response and excitation vectors between the two connection coordinates, respectively.The response and excitation of coordinatec1 follow this rule and are expressed as follows:
(60)Xc1II=HcFc1II,Hc=s2Mb−1,Mb=diagmbmbmbJbrJbrJbz,where mb is the mass of the bracket and Jbr and Jbz are the radial inertia and polar inertia, respectively.The transitive relation between response and excitation for the two connection coordinates is derived from (59) and (60) that is given as follows:
(61)Xc1IIXc2II=HcHcB1B1THcB1THcB1Fc1IIFc2II.Therefore, we obtain the FRF matrix of substructureII as follows:
(62)HII=HcHcB1B1THcB1THcB1.


## 4.4. Substructure III

Substructure III is the flexible interface that is demonstrated by an AHSP. Substructure III only considers a connection coordinate (coordinateb2) for the study and mainly investigates the dynamic responses of the CMG mounting point. The FRF matrix of substructure III is obtained by its finite element model (FEM).The structure of the AHSP consists of two layers of aluminum skins on the two sides with honeycomb core in the middle. The honeycomb core is equivalent to homogeneous orthotropic material, and thus the AHSP is considered a laminate with three-layer plates. The size of the AHSP is 1 m × 0.5 m × 0.03 m, and the thickness of a single-layer aluminum skin is 0.3 mm. Additionally,E, G, ρ, and υ denote Young’s modulus, shear modulus, density, and Poisson ratio, respectively. The material properties of the aluminum skins and the honeycomb core equivalent model are shown in Tables 1 and 2, respectively, and are derived from reference [29].Table 1
Material properties of the aluminum alloy.
Young’s modulus (E)
Density (ρ)
Poisson ratio (υ)
71 GPa
2700 kg/m3
0.33Table 2
Equivalent material properties of the honeycomb core.
Symbol
Values
Unit
Exx
0.0354
MPa
Eyy
0.0354
MPa
Ezz
655.87
MPa
Gxy
0.0266
MPa
Gyz
141.12
MPa
Gzx
92.463
MPa
ρ
24.94
kg/m3
υ12
0.999856
—
υ23
0
—
υ31
0
—Furthermore, MSC PATRAN is used to establish the FEM of the AHSP, as shown in Figure12. In the FEM, a laminate element is applied to mesh the AHSP. The equivalent material properties of the honeycomb core are simulated using a 2D orthotropic material. The nodes on the four sides of the FEM are imposed by six-DOF displacement constraints, and the separation distance between two nodes is 0.1 m. The FRF matrix of substructure III Hb2 is obtained via harmonic response analysis of the FEM by acting six-DOF unit force on the CMG mounting point independently.Figure 12
Schematic diagram of FEM of the AHSP.


## 4.5. FRF Matrix of the Coupling System

With respect to the substructure system, connection coordinates on the two sides of connection substructure 1 are not in agreement. In the process of providing control moments, the gimbal rotates corresponding to the bracket around the shaft, and thus coordinateg rotates to coordinate 0 around the y-axis as shown in Figure 13 in which θ is the rotation angle. Table 3 lists the angle relation between coordinates g and 0 in which the orientation cosine matrix is obtained as follows:
(63)A1=A300A3,A3=cosθ0cos0.5π−θ010cos0.5π+θ0cosθ.Figure 13
Angle relation between coordinatesg and 0.Table 3
Angle relation between the coordinatesg and 0.
Axis
xg
yg
zg
x0
θ
0.5π
0.5π − θ
y0
0.5π
0
0.5π
z0
0.5π + θ
0.5π
θConnection coordinates on the two sides of connection substructure 2 are accordant, and thus the orientation cosine matrix between coordinatesa and 0 is the unit matrix. The sum of the DOFs of all the internal and connection coordinates of the substructure system is six and twelve, respectively. Subsequently, the connection coordinate transformation matrix of the substructure system is obtained and expressed as follows:
(64)AF=diagI6A1I6I12.Finally, the FRF matrix of the substructure system is achieved by substituting (58), (62), and (64) into (19) and is expressed as follows:
(65)Moreover, (2), (33), and (34) are used, and the connection receptance matrix of connection substructure 1 is expressed as follows:
(66)R1o0x0y0z0=Z1−1=K1+sC1−1.Connection receptance of connection substructure 2 tends towards zero because it is a rigid connection. Hence, the connection receptance matrix of the substructure system is expressed as follows:(67)RF=Z1−1000.Finally, the FRF matrix of the coupling systemHF is achieved by substituting (65) and (67) into (15).


## 4.6. Dynamic Responses of the Coupling System

For the coupling system, only external excitation acts on the flywheel as discussed in Section4.2.4. and shown in (50). Thus, the external excitation vector of the coupling system is obtained as follows:
(68)fF=fsTfdT0TT.We apply the FT to (68), and the external excitation vector is transformed into the frequency domain form that is expressed as follows:
(69)FFs=FsTsFdTs0TT.Displacement response of the coupling system is obtained by using its FRF matrix to multiply the external excitation vector as follows:(70)XF=HFFF.The displacement response is multiplied bys2, and we obtain the acceleration response of the coupling system as follows:
(71)X¨F=HFFFs2.Reaction forces and moments between CMG and the flexible interface are achieved by the impedance matrix of the CMG mounting point multiplied by its displacement response vector as follows:(72)Fm=Hb1−1Xm,where Xm (a part of XF) is the displacement response vector of the CMG mounting point.


## 5. Dynamic Responses of CMG on the Rigid Interface

The FRF matrix of CMG on the rigid interface is obtained by the similar method with the flexible interface. As shown in Figure14, the substructure system of CMG on the rigid interface contains two substructures and a connection substructure in which definitions are essentially consistent with the flexible interface with the exception of the CMG mounting point, which is considered an internal coordinate (coordinate h) that exhibits zero displacement in this situation because CMG is hard-mounted on the foundation.Figure 14
Substructure division of CMG on the rigid interface.Based on the previous section, the FRF matrix of this substructure system is expressed as follows:(73)The connection receptance matrix of the substructure system isR1. Subsequently, the FRF matrix of CMG on the rigid interface HR is obtained by substituting (66) and (73) into (15).The reaction forces and moments between CMG and the rigid interface are considered external excitation acting on the CMG mounting point. Subsequently, the excitation vector of the system is expressed as follows:(74)FR=FfTFhT0TT,where Fh denotes the reaction forces and moments between CMG and the rigid interface.The FRF matrixHR is partitioned in internal coordinates f and h and connection coordinate c, and thus the relationship between the external excitation and displacement response vectors is expressed as follows:
(75)HRFR=XR,HffHfhHfcHhfHhhHhcHcfHchHccFfFh0=XfXhXc.With respect to the rigid interface, the internal coordinateh is restrained to zero displacement, and thus Xh=0, and the rigid interface reaction force and moment vector Fh is calculated as follows:
(76)Fh=−Hhh−1HhfFf.


## 6. Simulation and Analysis

In this section, numerical and multibody dynamics simulations are used to examine the dynamic characteristics of CMG on the flexible interface.


### 6.1. Numerical Simulation and Analysis

MATLAB numerical simulations are performed to simulate the reaction forces and moments of CMG on both flexible and rigid interfaces by using (72) and (76), respectively. Subsequently, the dynamic response characteristics of the coupling system are analyzed by comparing the simulation results of the two interfaces. The reaction forces and moments are simulated when the gimbal angle θ = l·30° (l=0,1,2,…,12) and the rotating speed of the flywheel stabilizes at 6000 rpm. The simulation parameters are shown in Table 4. Simulation frequency range and frequency precision are 0–300 Hz and 0.0625 Hz, respectively.Table 4
Simulation parameters of the coupling system.
Symbol
Description
Value
Unit
mf
Mass of the flywheel
8
kg
Jfr
Radial inertia of the flywheel
0.02
kg·m2
Jfz
Polar inertia of the flywheel
0.03
kg·m2
mg
Mass of the gimbal
10
kg
Jgr
Radial inertia of the gimbal
0.06
kg·m2
Jgz
Polar inertia of the gimbal
0.09
kg·m2
mb
Mass of the bracket
14
kg
Jbr
Radial inertia of the bracket
0.1
kg·m2
Jbz
Polar inertia of the bracket
0.2
kg·m2
kr0
Radial stiffness of the flywheel bearing
4 × 106
N/m
ka0
Axial stiffness of the flywheel bearing
2 × 107
N/m
cr0
Radial damping of the flywheel bearing
2000
N·s/m
ca0
Axial damping of the flywheel bearing
3000
N·s/m
kr1
Radial stiffness of the gimbal bearing
8 × 106
N/m
ka1
Axial stiffness of the gimbal bearing
4 × 107
N/m
cr1
Radial damping of the gimbal bearing
2000
N·s/m
ca1
Axial damping of the gimbal bearing
3000
N·s/m
ks
Servo dynamic stiffness of the servo system
8 × 106
N·m/rad
d0
Distance from origino1 to the flywheel bearing
0.04
m
d1
Distance from origino2 to the gimbal bearing
0.055
m
h
Distance between connection coordinatesc1 and o1
0.06
m
Us
Static mass imbalances
3.6 × 10−6
kg·m
Ud
Dynamic mass imbalances
7.02 × 10−8
kg·m2
φs
Initial phase of the static mass imbalance
0
rad
φd
Initial phase of the dynamic mass imbalance
0
rad
Ω
Rotating speed of the flywheel
200π
rad/sFigure15 provides the simulation results of the force and moment amplitude of the two interfaces at the fundamental frequency 100 Hz with respect to variations in the gimbal angle. Figure 15 indicates that the dynamic characteristics of the flexible interface are considerably different from those of the rigid interface. With respect to the six-component amplitude of the rigid interface, Fy and Mx are essentially invariable, and the extreme values of the other components are located in θ = l·90° (l=1,2,…,5) and their change rules approximate the sine curve. Conversely, those of the flexible interface did not follow these rules in which all the components vary with the gimbal angle, and the extreme values of Fz shifts 30° to the right of the abscissa and Mz presents the reverse trend at 90° and 270°. Additionally, disturbance amplitudes of the flexible interface are evidently higher than those of the rigid interface with the exception of Fz, and the mean values of these disturbance amplitudes on both interfaces and their ratio are shown in Table 5. The table indicates that the mean values of the three force components of the rigid interface are extremely close to each other while those of the flexible interface exhibit a significant difference, Fz, which is considerably lower.Figure 15
Force and moment amplitude on both the flexible interface and rigid interface.Table 5
Mean values of disturbance amplitude on two interfaces.
Component
Flexible interface
Rigid interface
Ratio
Fx
3.99 N
2.89 N
1.38
Fy
5.84 N
2.45 N
2.38
Fz
0.64 N
2.54 N
0.25
Mx
0.57 N·m
0.15 N·m
3.81
My
0.31 N·m
0.18 N·m
1.72
Mz
0.03 N·m
0.01 N·m
2.11


### 6.2. Multibody Dynamics Simulation Validation

The dynamic responses of the coupling system on the DSM are validated by performing both numerical and multibody dynamics simulations to simulate acceleration response of the CMG mounting point. First, MATLAB simulates the acceleration response via (64). Subsequently, MSC Adams conducts multibody dynamics simulation of the coupling system. Finally, the simulation results of the two simulations are compared, and the dynamic responses of the coupling system on the DSM are validated.In the Adams dynamic model, the flywheel, gimbal, and bracket are simulated by the rigid body parts, the bearings by the flexibility characteristic bearing parts, and the servo dynamic stiffness by the torsion spring part. The flexible interface is established by importing a modal neutral file (MNF) of the AHSP that is exported from modal analysis of the AHSP in PATRAN. The imbalance mass is attached to the flywheel, and the rotating speed drive is applied to perform the Adams model around thezg-axis. The sampling frequency and sampling time of the Adams simulations are 2048 Hz and 16 s, respectively. The simulation results of the time domain are processed and transformed into the frequency domain via FT.Simulation results of the translational acceleration (TA) and angular acceleration (AA) amplitude of the two simulations at the fundamental frequency of 100 Hz vary with the gimbal as shown in Figure16. The figure indicates that the two simulation results are essentially consistent. The simulation error of the dynamic responses on the DSM is within ±6% as shown in Table 6. The results in Figure 16 and Table 6 indicate that the DSM accurately predicts the coupled microvibrations of CMG on the flexible interface.Figure 16
Translational and angular acceleration amplitude by MATLAB and Adams simulations.Table 6
Simulation error of the dynamic response on the DSM (%).
Gimbal angle
TAx
TAy
TAz
AAx
AAy
AAz
0°
−4.44
−3.47
−3.31
−3.15
−4.48
−3.47
30°
−3.96
−3.63
−2.23
−3.47
−3.84
−3.91
60°
−2.96
−4.08
−1.25
−4.21
−2.64
−4.51
90°
−2.11
−4.37
−0.81
−4.74
−1.08
−4.69
120°
−2.75
−3.84
−1.12
−4.20
−4.23
−3.80
150°
−3.73
−3.21
−2.17
−3.45
−4.69
−3.10
180°
−4.50
−2.93
−4.01
−3.02
−5.15
−2.94
210°
−4.01
−3.19
−2.09
−3.42
−4.35
−3.30
240°
−2.91
−3.88
−1.20
−4.31
−2.79
−3.82
270°
−1.87
−4.37
−0.80
−4.76
−0.83
−3.74
300°
−2.84
−4.03
−1.12
−4.06
−4.09
−2.98
330°
−3.69
−3.64
−2.11
−3.40
−3.90
−3.16
360°
−4.44
−3.47
−3.31
−3.16
−4.49
−3.47


## 6.1. Numerical Simulation and Analysis

MATLAB numerical simulations are performed to simulate the reaction forces and moments of CMG on both flexible and rigid interfaces by using (72) and (76), respectively. Subsequently, the dynamic response characteristics of the coupling system are analyzed by comparing the simulation results of the two interfaces. The reaction forces and moments are simulated when the gimbal angle θ = l·30° (l=0,1,2,…,12) and the rotating speed of the flywheel stabilizes at 6000 rpm. The simulation parameters are shown in Table 4. Simulation frequency range and frequency precision are 0–300 Hz and 0.0625 Hz, respectively.Table 4
Simulation parameters of the coupling system.
Symbol
Description
Value
Unit
mf
Mass of the flywheel
8
kg
Jfr
Radial inertia of the flywheel
0.02
kg·m2
Jfz
Polar inertia of the flywheel
0.03
kg·m2
mg
Mass of the gimbal
10
kg
Jgr
Radial inertia of the gimbal
0.06
kg·m2
Jgz
Polar inertia of the gimbal
0.09
kg·m2
mb
Mass of the bracket
14
kg
Jbr
Radial inertia of the bracket
0.1
kg·m2
Jbz
Polar inertia of the bracket
0.2
kg·m2
kr0
Radial stiffness of the flywheel bearing
4 × 106
N/m
ka0
Axial stiffness of the flywheel bearing
2 × 107
N/m
cr0
Radial damping of the flywheel bearing
2000
N·s/m
ca0
Axial damping of the flywheel bearing
3000
N·s/m
kr1
Radial stiffness of the gimbal bearing
8 × 106
N/m
ka1
Axial stiffness of the gimbal bearing
4 × 107
N/m
cr1
Radial damping of the gimbal bearing
2000
N·s/m
ca1
Axial damping of the gimbal bearing
3000
N·s/m
ks
Servo dynamic stiffness of the servo system
8 × 106
N·m/rad
d0
Distance from origino1 to the flywheel bearing
0.04
m
d1
Distance from origino2 to the gimbal bearing
0.055
m
h
Distance between connection coordinatesc1 and o1
0.06
m
Us
Static mass imbalances
3.6 × 10−6
kg·m
Ud
Dynamic mass imbalances
7.02 × 10−8
kg·m2
φs
Initial phase of the static mass imbalance
0
rad
φd
Initial phase of the dynamic mass imbalance
0
rad
Ω
Rotating speed of the flywheel
200π
rad/sFigure15 provides the simulation results of the force and moment amplitude of the two interfaces at the fundamental frequency 100 Hz with respect to variations in the gimbal angle. Figure 15 indicates that the dynamic characteristics of the flexible interface are considerably different from those of the rigid interface. With respect to the six-component amplitude of the rigid interface, Fy and Mx are essentially invariable, and the extreme values of the other components are located in θ = l·90° (l=1,2,…,5) and their change rules approximate the sine curve. Conversely, those of the flexible interface did not follow these rules in which all the components vary with the gimbal angle, and the extreme values of Fz shifts 30° to the right of the abscissa and Mz presents the reverse trend at 90° and 270°. Additionally, disturbance amplitudes of the flexible interface are evidently higher than those of the rigid interface with the exception of Fz, and the mean values of these disturbance amplitudes on both interfaces and their ratio are shown in Table 5. The table indicates that the mean values of the three force components of the rigid interface are extremely close to each other while those of the flexible interface exhibit a significant difference, Fz, which is considerably lower.Figure 15
Force and moment amplitude on both the flexible interface and rigid interface.Table 5
Mean values of disturbance amplitude on two interfaces.
Component
Flexible interface
Rigid interface
Ratio
Fx
3.99 N
2.89 N
1.38
Fy
5.84 N
2.45 N
2.38
Fz
0.64 N
2.54 N
0.25
Mx
0.57 N·m
0.15 N·m
3.81
My
0.31 N·m
0.18 N·m
1.72
Mz
0.03 N·m
0.01 N·m
2.11


## 6.2. Multibody Dynamics Simulation Validation

The dynamic responses of the coupling system on the DSM are validated by performing both numerical and multibody dynamics simulations to simulate acceleration response of the CMG mounting point. First, MATLAB simulates the acceleration response via (64). Subsequently, MSC Adams conducts multibody dynamics simulation of the coupling system. Finally, the simulation results of the two simulations are compared, and the dynamic responses of the coupling system on the DSM are validated.In the Adams dynamic model, the flywheel, gimbal, and bracket are simulated by the rigid body parts, the bearings by the flexibility characteristic bearing parts, and the servo dynamic stiffness by the torsion spring part. The flexible interface is established by importing a modal neutral file (MNF) of the AHSP that is exported from modal analysis of the AHSP in PATRAN. The imbalance mass is attached to the flywheel, and the rotating speed drive is applied to perform the Adams model around thezg-axis. The sampling frequency and sampling time of the Adams simulations are 2048 Hz and 16 s, respectively. The simulation results of the time domain are processed and transformed into the frequency domain via FT.Simulation results of the translational acceleration (TA) and angular acceleration (AA) amplitude of the two simulations at the fundamental frequency of 100 Hz vary with the gimbal as shown in Figure16. The figure indicates that the two simulation results are essentially consistent. The simulation error of the dynamic responses on the DSM is within ±6% as shown in Table 6. The results in Figure 16 and Table 6 indicate that the DSM accurately predicts the coupled microvibrations of CMG on the flexible interface.Figure 16
Translational and angular acceleration amplitude by MATLAB and Adams simulations.Table 6
Simulation error of the dynamic response on the DSM (%).
Gimbal angle
TAx
TAy
TAz
AAx
AAy
AAz
0°
−4.44
−3.47
−3.31
−3.15
−4.48
−3.47
30°
−3.96
−3.63
−2.23
−3.47
−3.84
−3.91
60°
−2.96
−4.08
−1.25
−4.21
−2.64
−4.51
90°
−2.11
−4.37
−0.81
−4.74
−1.08
−4.69
120°
−2.75
−3.84
−1.12
−4.20
−4.23
−3.80
150°
−3.73
−3.21
−2.17
−3.45
−4.69
−3.10
180°
−4.50
−2.93
−4.01
−3.02
−5.15
−2.94
210°
−4.01
−3.19
−2.09
−3.42
−4.35
−3.30
240°
−2.91
−3.88
−1.20
−4.31
−2.79
−3.82
270°
−1.87
−4.37
−0.80
−4.76
−0.83
−3.74
300°
−2.84
−4.03
−1.12
−4.06
−4.09
−2.98
330°
−3.69
−3.64
−2.11
−3.40
−3.90
−3.16
360°
−4.44
−3.47
−3.31
−3.16
−4.49
−3.47


## 7. Conclusion

The aim of the study involves investigating the microvibration characteristics of the CMG-flexible interface coupling system by the DSM. Initially, a DSM based on FRF is established that is extremely suitable for simultaneously tackling multiple different coordinate substructures with hybrid connections and exhibits high efficiency. Subsequently, the method is employed to obtain the FRF matrix of CMG on both flexible and rigid interfaces, and this is used to predict the reaction forces and moments between CMG and the mounting interface and the acceleration response of the CMG mounting point. MATLAB simulations are conducted to simulate the reaction forces and moments on the two interfaces. The simulation results indicate that evident differences exist between the dynamic characteristics of the two interfaces regardless of the locations of the extreme value or the change rules of the amplitude. Additionally, the amplitudes of the flexible interface evidently exceed those of the rigid interface with the exception ofFz. The results indicate that the traditional coupling disturbance analysis method that applies the CMG hard-mount testing spectrum to a spacecraft FRF to predict its performance produces a high error. Finally, multibody dynamics simulations are performed to validate the dynamic responses of the coupling system achieved by the established DSM and suggest results wherein predictions of the FRF matrix achieved by the DSM consistently predicted the results of multibody dynamics simulations with a simulation error of less than six percent. Therefore, the DSM is used to predict coupled microvibration characteristics between the CMG and flexible interface and exhibits high prediction accuracy and computational efficiency.



---

*Source: 5045740-2018-04-08.xml*