**Analytic Model and the Influence of Actuator Number on the Performance of Plasma Synthetic Jet Actuator Array**

**Shengfang Huang 1, Zhibo Zhang 1,\*, Huimin Song 1,\*, Yun Wu 2, Zhengzhong Sun <sup>3</sup> and Yinghong Li <sup>1</sup>**


Received: 6 August 2018; Accepted: 29 August 2018; Published: 3 September 2018

**Abstract:** Coupled with the multichannel discharge model and plasma synthetic jet actuator (PSJA) aerodynamic model, an analytical model to predict the performance of the PSJA array is put forward. The multichannel discharge model takes these factors into consideration, the delay time in the breakdown process, the electrical transformation of the discharge channel from a capacitor to a resistor induced by the air breakdown, and the varying plasma resistance in the discharge process. The PSJA aerodynamic model is developed based on the conservation equations of mass, momentum, energy, and the lumped capacitance method. The multichannel discharge model can simulate the multichannel discharge process and give the discharge energy in the plasma channel. With a constant heating efficiency, the time-independent heating energy deposition power in a discharge channel is obtained. Importing the heating energy, the PSJA aerodynamic model presents the evolution process of the jet. Simulation results show that the jet strength induced by a single actuator decreases with the number of actuators in the PSJA array. When the actuator number increases from 1 to 20, the weakening extent of mass ejected, peak jet velocity, and jet duration time is 62%, 54%, and 33%, respectively. The discharge efficiency increases with the actuator number, while the thermodynamic efficiency decreases with the actuator number. As a result, the total energy efficiency doesn't always increase with an increase in the number of actuators. When the discharge efficiency of a conventional one channel discharge has been a relatively large value, the total energy efficiency actually decreases with the growth of actuator number.

**Keywords:** Plasma flow control; multichannel discharge; plasma synthetic actuator; actuator array; analytic model

## **1. Introduction**

Plasma actuators have attracted a lot of research attention in recent years because of the associated advantages, such as the absence of moving components, fast response, and wide bandwidth. So far promising control outcomes have been revealed in suppressing flow separation [1,2], alleviating turbomachinery stall [3,4], and shock wave/boundary layer interaction [5,6]. A plasma synthetic jet actuator (PSJA) is one of the typical plasma actuators. The PSJA makes use of the heat released from the spark discharge process between electrodes within a small cavity where the electrodes are

housed. As the velocity of the resulting jet is high (larger than 100 m/s), the PSJA is suitable for high-speed applications. It should be noted that the dielectric barrier discharge (DBD) plasma actuator with annular electrodes can also generate a vertical jet [7,8], and has the name of plasma synthetic jet actuator too [9]. However, the jet velocity from DBD-PSJA is much less. The present work has focused on the PSJA with the principle of spark discharge.

As the diameter of the PSJA orifice is in the order of millimeters, to enlarge the affected area, an array containing distributed actuators is necessary for practical applications [10]. In 2009, Caruana et al. [11] attempted to reduce the jet exhaust noise using six PSJAs. Unfortunately, little effect was achieved, and it was concluded that this was as a result of an insufficient number of actuators. In 2012, an array of three PSJAs was used to control the unsteadiness of a shock wave/turbulent boundary layer by Narayanaswamy et al. [6]. In 2013, an attempt of using 20 PSJAs is made to reduce the trailing edge separation on a NACA-0015 profile [12,13]. Therefore, it can be found that there is a strong desire to implement more PSJAs in various flow controls to deliver a strong actuation.

The arc discharge channel was revealed to feature negative resistant characteristics [14,15], namely, an increase in the current leads to a decrease in the voltage. As a result, the multichannel discharge circuit is difficult to power when using a single power supply. In 2005, M. Samimy et al. [16] tested the discharge in eight channels connected in parallel by using a 30 kΩ current-limiting resistor. This large current-limiting resistor results in a large energy waste. As introduced above, 20 PSJAs were driven simultaneously, but 20 high power supply systems are used by Caruana [12], which is apparently not suitable for on-board applications. Recently, a multichannel discharge technique based on a voltage relay concept is proposed by the authors' group [17–19]. It allows multi-channel discharge without increasing the input voltage and without adding an additional current-limiting resistor. The new voltage relay circuit is a breakthrough. Despite its initial success in experimental testing, the manner in which the channel number affects the jet flow is still not clear, which is of great importance in optimizing the PSJA array's aerodynamic performance.

One efficient way to study the effect of actuator number is to use an analytical model. In 2003, Grossman et al. [20] developed a first-order PSJA model. However, it can only predict the basic performance parameter, such as the jet velocity. In 2010, Haack et al. [21] made a further development where the jet generation process is included, but without modeling the recovery process. In 2015, Zong et al. [22] developed a model that can predict a PSJA's entire working cycle. However, the models introduced above only simulate the operation of a single PSJA. No model is available to predict the performance of a PSJA array, namely a number of PSJAs.

The plasma synthetic jet actuator array (PSJA array) developed by the authors can enlarge the area affected by the actuator greatly, which is more suitable for practical application. A new analytic model is developed in the present work. It consists of two parts: The multichannel discharge model and the PSJA aerodynamic model. The new model enables simulation of the whole working cycle of the PSJA array, and allows the comparison with that of a single PSJA. As the number of PSJAs potentially affects the performance, the effect of the number of PSJAs is also studied using the new model.

#### **2. Analytic Model Development**

Same as the conventional PSJA, the working process of the PSJA array consists of three stages: Energy deposition, jet generation, and recovery. However, owing to the disparity in the discharge circuit, the energy deposition process is different. After the energy deposition process, the jet generation and recovery process are similar, which can be described by conservation of mass, momentum, and energy equations. Hence, the analytical model of the PSJA array is divided into two models: A multichannel discharge model and PSJA aerodynamic model. The multichannel discharge model is used to describe the discharge process and calculate the energy deposited in each discharge channel. The PSJA model is used to describe the jet generation and recovery process. The deposited energy of PSJA is obtained from the multichannel discharge model. The multichannel discharge model has been described in reference [23].

## *2.1. The PSJA Aerodynamic Model*

As shown in Figure 1, the PSJA consists of four main parts: Throat, chamber, shell, and electrodes. To simplify the PSJA aerodynamic model, the fluid region including the throat and chamber is treated as two control volumes. The pressure, temperature, and density are averaged over the control volume. Then, the conservation of mass, momentum, and energy equations are manipulated to describe the changes in density, jet velocity, pressure, and temperature in the control volume.

**Figure 1.** The structure of plasma synthetic jet actuator (PSJA).

To simplify these conservation equations, some assumptions are adopted in this model.


## *2.2. The Conservation Equations*

For inviscid flow, the integral form of the conservation equations is given as following:

$$\frac{\partial}{\partial t} \iiint\_{v} \rho \, d\upsilon + \iint\_{s} \rho \stackrel{\rightarrow}{\boldsymbol{\mu}} \stackrel{\rightarrow}{d} \stackrel{\rightarrow}{S} = 0 \tag{1}$$

$$\frac{\partial}{\partial t} \iiint\_{\upsilon} \rho \stackrel{\rightarrow}{u} dv + \iint\_{s} \left( \rho \stackrel{\rightarrow}{u} \stackrel{\rightarrow}{d} \stackrel{\rightarrow}{S} \right) \stackrel{\rightarrow}{u} = - \iint\_{s} P \stackrel{\rightarrow}{S} + \iiint\_{\upsilon} \rho \stackrel{\rightarrow}{f} dv \tag{2}$$

$$\frac{\partial}{\partial t} \iiint\_{\upsilon} \rho \left(\varepsilon + \frac{\mu^2}{2}\right) d\upsilon + \iint\_{\varsigma} \rho \left(\varepsilon + \frac{\mu^2}{2}\right) \overrightarrow{\bar{u}} \, d\bar{S} = -\iint\_{\varsigma} P \overrightarrow{\bar{u}} \, d\bar{S} + \dot{Q} \tag{3}$$

For the chamber control volume, the gas velocity is negligible. These conservation equations are represented as the following:

$$\frac{d\rho\_{\rm c}}{dt} = -\frac{\rho\_{\rm i}u\_{\rm i}A\_{\rm i}}{V\_{\rm c}}\tag{4}$$

$$\frac{dT\_{\rm c}}{dt} = \frac{1}{\rho\_{\rm c}} \left( \frac{\dot{Q} - C\_{\rm v} T\_{\rm i} \rho\_{\rm i} \mu A\_{\rm t} - 0.5 \rho\_{\rm i} \mu^3 A\_{\rm t}}{C\_{\rm v} V\_{\rm c}} - T\_{\rm c} \frac{d \rho\_{\rm c}}{dt} \right) \tag{5}$$

For the throat control volume, these conservation equations are represented as the following:

$$\frac{d\rho\_{\rm t}}{dt} = \frac{\rho\_{\rm i}u\_{\rm i}A\_{\rm i} - \rho\_{\rm o}u\_{\rm o}A\_{\rm o}}{V\_{\rm t}}\tag{6}$$

$$\frac{\mathbf{d}u}{\mathbf{d}t} = \rho\_\mathbf{t} \left( \frac{P\_\mathbf{i}A\_\mathbf{i} - P\_\mathbf{o}A\_\mathbf{o} + \rho\_\mathbf{i}u^2A\_\mathbf{i} - \rho\_\mathbf{o}u^2\_\mathbf{o}A\_\mathbf{o}}{V\_\mathbf{t}} - u\frac{\mathbf{d}\rho\_\mathbf{t}}{\mathbf{d}t} \right) \tag{7}$$

$$\frac{dT\_{\rm l}}{dt} = \frac{1}{\rho\_{\rm l}} \left( \frac{(P\_{\rm l} - P\_{\rm o})\mu A\_{\rm l} + \left( \left( \mathbb{C} \gamma (T\_{\rm l} \rho\_{\rm l} - T\_{\rm o} \rho\_{\rm o}) + 0.5 (\rho\_{\rm l} - \rho\_{\rm o})u^{2} \right) \mu A\_{\rm l} \right) - V\_{\rm l} \left( 2\rho\_{\rm l} u \frac{du}{dt} + u^{2} \frac{d\rho\_{\rm l}}{dt} \right)}{\mathbb{C}\_{\rm v} V\_{\rm l}} - T\_{\rm l} \frac{d\rho\_{\rm l}}{dt} \right) \tag{8}$$

In the Equations (4)–(8), the subscripts c and t represent the averaged parameters over the chamber control volume and the throat control volume, respectively; the subscripts i and o represent the averaged parameters over the throat inlet control face and the throat outlet control face, respectively. The *u* indicates the averaged gas velocity over the throat control volume.

To solve these equations, these parameters including the pressure, density, and temperature over the throat inlet and outlet control face must be calculated. As shown in Figure 2, in a different working stage, the gas flow direction is the opposite. As a result, the way to calculate the parameters over the throat inlet and outlet control face is different.

**Figure 2.** The schematic diagram of the two working stages of PSJA.

In jet stage, the gas flow from the chamber to the outside environment. At this time, the parameters over the throat inlet control face are determined by the gas state parameters in the chamber, and the parameters over the throat outlet control face are determined by the gas state parameters in the throat. Additionally, when the gas flows through the throat, the diameter of the pipeline decreases quickly, which induced the local pressure loss. In this model, the local pressure loss is considered using an engineering method, as shown in Equation (9).

$$P\_{\rm loss} = 0.5 \rho u^2 \ast 0.5 \left( 1 - \frac{A\_2}{A\_1} \right) + \Delta p\_{pl} \tag{9}$$

Taking these factors into consideration, the corresponding functions are presented as the following:

$$\begin{cases} \rho\_{\text{i}} = \rho\_{\text{c}} \varepsilon(\lambda) \\ P\_{\text{i}} = P\_{\text{c}} \pi(\lambda) - 0.5 \rho\_{\text{i}} u^{2} 0.5 \left( 1 - \frac{A\_{\text{i}}}{A\_{\text{c}}} \right) - \Delta p\_{pl} \\ T\_{\text{i}} = \frac{P\_{\text{i}}}{R \rho\_{\text{i}}} \\ \rho\_{\text{o}} = \rho\_{\text{t}} \\ P\_{\text{o}} = P\_{\text{os}} \\ T\_{\text{o}} = \frac{P\_{\text{o}}}{R \rho\_{\text{o}}} \end{cases} \tag{10}$$

*Appl. Sci.* **2018**, *8*, 1534

However, limiting by the throat geometry, the maximum gas velocity is the local sound speed. That means that the flow condition has two modes: To be chocked and to be unchocked. The criterion for judgment is as follows:

$$\begin{cases} \begin{array}{c} P\_{\text{i}} + 0.5 \rho\_{\text{i}} u^{2} \geq 1.89 P\_{\infty} \\ \gamma \geq 1 \end{array} \tag{11} $$

If flow parameters meet this criterion, the flow is considered to be choked. The pressure at the throat outlet face is calculated as follows:

$$P\_0 = \frac{1}{1.89} \left( P\_{\text{l}} + 0.5 \rho\_{\text{l}} u^2 \right) \tag{12}$$

In recovery stage, the gas flows from the outside environment to the chamber. At this time, the parameters over the throat inlet control face are determined by the gas state parameters in the throat; the parameters over the throat outlet control face are determined by the gas state parameters in the external environment. The corresponding functions are presented as follows:

$$\begin{cases} \rho\_{\text{i}} = \rho\_{\text{t}}\\ P\_{\text{i}} = P\_{\text{c}}\\ T\_{\text{i}} = \frac{P\_{\text{i}}}{R\rho\_{\text{i}}}\\ \rho\_{\text{o}} = \rho\_{\text{os}}\varepsilon(\lambda)\\ P\_{\text{i}} = P\_{\text{os}}\pi(\lambda) - 0.25\rho\_{\text{o}}\mu^{2} \\ T\_{\text{i}} = \frac{P\_{\text{o}}}{R\rho\_{\text{o}}} \end{cases} \tag{13}$$

In the above functions, *λ* indicates the velocity coefficient, which is calculated by the following functions:

$$\begin{cases} a\_{\rm cr} = \sqrt{\frac{2\gamma}{\gamma+1}RT\_t^\*} \\ T\_t^\* = T\_t + \frac{\mu^2}{2\Gamma\_p} \\ \lambda = \frac{|\mu|}{a\_{\rm cr}} \end{cases} \tag{14}$$

*ε*(*λ*) and *π*(*λ*) are calculated as follows:

$$\begin{aligned} \varepsilon(\lambda) &= \left(1 - \frac{\gamma - 1}{\gamma + 1} \lambda^2\right)^{\frac{1}{\gamma - 1}} \\ \pi(\lambda) &= \left(1 - \frac{\gamma - 1}{\gamma + 1} \lambda^2\right)^{\frac{\gamma}{\gamma - 1}} \end{aligned} \tag{15}$$

## *2.3. The Thermal Modeling*

In Equation (5), *Q* represents the discharge heating power and the thermal energy loss. The discharge heating power can be represented by Ohm heating easily. The thermal energy loss through the wall and the electrode is modeled based on the lumped capacitance method.

The heat transfer process in the PSJA is shown in Figure 3. Firstly, by forced convection, the heat transfers from the chamber gas to the shell inner wall. Secondly, by heat conduction through the electrodes and shell wall, the heat transfers from the shell inner wall to the shell outer wall. Thirdly, by natural convection, the heat transfers from the shell outer wall to the external air.

Similar to the electric circuit, the heat transfer process can be described as a thermal circuit, as shown in Figure 4. *T*c, *T*w,in, *T*w,out, and *T*<sup>∞</sup> indicate the temperature of chamber gas, shell inner wall, shell outer wall, and external air, respectively. *R*in, *R*e, *R*w, and *R*out represent the thermal resistance of force convection, electrode conduction, shell conduction, and natural convection. *C*<sup>w</sup> represents the thermal capacitance of the shell.

**Figure 3.** The schematic diagram of the heat transfer process.

**Figure 4.** Electrical representation of the thermal heat transfer process.

Based on the theory of engineering thermal transmission, the thermal resistance of force convection, electrode conduction, and natural convection is given as follows.

$$\begin{cases} \begin{aligned} R\_{\text{in}} &= \frac{1}{l\_{\text{in}} A\_{\text{c,in}}}\\ R\_{\text{out}} &= \frac{1}{l\_{\text{out}} A\_{\text{c,out}}}\\ R\_{\text{e}} &= \frac{0.5 l\_{\text{e}}}{\kappa\_{\text{e}} A\_{\text{e}}} \end{aligned} \end{cases} \tag{16}$$

The structure of the shell wall is complex, and is regarded as a combination of three simple structures, as shown in Figure 5. The upper structure of the shell is seen as a single wall. The middle structure of the shell is seen as a single hollow cylinder wall. The subjacent structure of the shell is also seen as a single wall. The thermal resistance of the three structures is given as follows.

$$\begin{cases} R\_{\text{up}} = \frac{h\_3}{\kappa\_{\text{up}} \frac{\pi}{4} \left( d\_1^2 - d\_3^2 \right)}\\ R\_{\text{out}} = R\_{\text{mid}} = \frac{\ln \frac{d\_2}{d\_1}}{2 \pi \kappa\_{\text{mid}} h\_1} \\ R\_{\text{sub}} = \frac{h\_2}{\kappa\_{\text{sub}} \frac{\pi}{4} d\_1^2} \end{cases} \tag{17}$$

The total thermal resistance of the shell wall is calculated as Equation (18).

$$R\_{\rm w} = \frac{1}{1/R\_{\rm up} + 1/R\_{\rm mid} + 1/R\_{\rm sub}}\tag{18}$$

With all of the circuit components defined, the differential function obtained from the circuit shown in Figure 4 is used to solve for the wall temperature.

$$\frac{dT\_{\rm w,in}}{dt} = \frac{1}{C\_{\rm w}} \left( \frac{T\_{\rm c} - T\_{\rm w,in}}{R\_{\rm in}} - \frac{T\_{\rm w,in} - T\_{\rm w}}{R\_{\rm out} + R\_{\rm c} R\_{\rm w} / \left(R\_{\rm c} + R\_{\rm w}\right)} \right) \tag{19}$$

Hence, combined with the discharge heating source, the parameter *Q* in Equation (10) is given.

$$\dot{Q} = Q\_{\rm arc} - \frac{T\_{\rm c} - T\_{\rm w,in}}{R\_{\rm in}} = 0.5 I(t)^2 R(t) - \frac{T\_{\rm c} - T\_{\rm w,in}}{R\_{\rm in}} \tag{20}$$

It must be pointed out that the energy calculated directly using Ohm heating theory is not the heating energy. The sheath energy loss and radiation energy loss must be taken into consideration. Based on the analysis of Guillaume Dufour [19], a heating efficiency of 50% is adopted.

**Figure 5.** The disassembly structure of the shell wall.

## *2.4. Model Validation*

To validate the model's reliability and accuracy, simulation results with ANSYS CFX software are compared with the results calculated from the PSJA aerodynamic model. The computation mesh created by ANSYS ICEM is fully structured, and the boundary type is labeled in Figure 6. The actuator geometry parameters are shown in Figure 7.

**Figure 6.** Computation grid.

**Figure 7.** The geometry parameters of actuator.

The solver is coupled and fully implicit with 2nd order accuracy in both time and space. A variable time step is used, which begins with 5 ns and is increased gradually to 1 μs in the

simulation. Considering the high turbulent vortex produced by PSJA, the SST turbulent model is chosen. The heating power is coupled into the Navier-Stokes equations as a heat source and expressed by a CFX expression. As the purpose of this step is only to validate the PSJA model, the discharge heating source is simplified as a constant value, as shown in Equation (21). *Q*<sup>h</sup> is the total heating energy; 0.005 J is adopted. *T* is the heating time; 100 ns is chosen. *V* is the volume of the chamber; approximately 88e−<sup>9</sup> m−3. Additionally, it is assumed the actuator is uniformly heated in this simulation.

$$\dot{Q} = \begin{cases} \frac{Q\_{\text{h}}}{TV} \text{ Wm}^{-3} \, t \le T\\ 0 \text{ Wm}^{-3} \, t \ge T \end{cases} \tag{21}$$

The jet velocity at the throat outlet calculated by the ANSYS CFX and PSJA aerodynamic model is shown in Figure 8a. The results obtained from the CFX with different grid number agree well. Therefore, the grid convergence has been reached in CFX. Moreover, the small difference in the results using different turbulent models suggests that the results are not sensitive to the turbulent model. The value obtained by ANSYS CFX is 141 m/s, and the value obtained by the PSJA model is 143 m/s. The difference is only 2 m/s. The disparity of jet termination time is also little, less than 3%. The value obtained by ANSYS CFX is 269 μs, and the value obtained by the PSJA model is 275 μs. The history of gauge pressure in the chamber obtained by the two methods is also similar, as shown in Figure 8b. Overall, the PSJA model can simulate the actuator characteristic with great accuracy.

**Figure 8.** The comparison between the CFD results and model results. (**a**) is the jet velocity at the throat outlet and (**b**) is the gauge pressure in the chamber.

## **3. Results and Discussion**

## *3.1. The Discharge Energy and Discharge Efficiency*

Based on the multichannel discharge technique, the discharge channel number can be increased greatly. When the energy stored in the discharge capacitor is fixed, the energy deposited in each channel varies with the discharge channel number. As energy is a critical factor determining the performance of the PSJA, the energy characteristic is first investigated.

The discharge capacitor is set as 10 nF, and the initial voltage is set as 6 kV. The discharge voltage of each electrode couple is set as 4 kV. The inductance is 1.65 μH. The sum of wire resistance and the equivalent serial resistance of the capacitor are chosen as 1.89 Ω. The discharge voltage and current waveforms with channel numbers of 1, 5, 10, and 20 are simulated. Based on the current and resistance of each discharge channel, the discharge power and discharge energy can be obtained.

The discharge voltage and current waveforms with different actuator numbers are plotted in Figure 9. As described in [16], the breakdown of each electrode couple happens in sequence, the time to form a complete discharge channel increases with the growth of the actuator number. What's more, with the increase in the actuator number, the total length of the discharge channel is enlarged, leading

to the growth of plasma resistance. As a result, the discharge current decreases. The deposited energy in each actuator ought to decrease.

**Figure 9.** The discharge waveforms with different actuator numbers. (**a**) shows the the voltage waveforms and (**b**) shows the current waveforms.

The energy deposition power and the discharge energy in the first PSJA are obtained from the simulation results, as shown in Figure 10. The energy deposition power (*P*deposition), and discharge energy (*Q*discharge) are defined in Equation (22). With the increase in the actuator number, the energy deposition power decreases. However, as shown in Figure 10b, the deposited energy doesn't decrease linearly. Since the gap distance between electrodes in all actuators is the same, and the actuator works in series, the resistance of the plasma channel in each actuator can be considered as the same. The discharge energy in each actuator can be expressed as Equation (23), where *Q* is the energy stored in the capacitor, *R*arc indicates the resistance of plasma channel, n is the actuator number, and *R*<sup>0</sup> is the external resistance in the discharge circuit. Therefore, if the actuator number increases, the deposited energy in one actuator doesn't decrease linearly. When the actuator number increases from 1 to 20, the deposited energy in the first actuator decreases 76.7% instead of 95%. In this way, the total energy deposited in all actuators increases, as shown in Figure 11. This means more energy stored in the capacitor is released in the discharge channel, leading to the improvement of discharge efficiency *η*d, which is defined in Equation (24).

$$Q\_{\text{discharge}} = \int P\_{\text{deposition}}(t)dt = \int I^2(t) / \lg(t)dt\tag{22}$$

$$Q\_{\text{discharge}} = Q \frac{R\_{\text{arc}}}{nR\_{\text{arc}} + R\_0} = \frac{Q}{n + R\_0/R\_{\text{arc}}} \tag{23}$$

$$\eta\_{\rm d} = \frac{\sum\_{k=1}^{k=n} Q\_{\rm discharge,k}}{Q} = \frac{\sum\_{k=1}^{k=n} \int I^2(t) / \mathcal{g}\_k(t) dt}{Q} \tag{24}$$

Figure 12 shows the deposited energy in each actuator with different actuator numbers. This figure shows that in this multichannel discharge mode, the deposited energy in a different actuator is almost the same. The coefficient of variation with 5, 10, and 20 actuators is 3.8%, 3.32%, and 4.77%, respectively. Based on the previous research [18], the breakdown of the electrode couple in each actuator happens in sequence. Coupled with the results shown in Figure 12, it concludes that although the time when the breakdown of the electrode couple in each actuator happens is different, the deposited energy is the same. The breakdown order has little influence on the deposited energy.

**Figure 10.** The discharge characteristics with different actuator numbers. (**a**) is the energy deposition power and (**b**) is the discharge energy in one actuator.

**Figure 11.** The total discharge energy and discharge efficiency versus actuator number.

**Figure 12.** The deposited energy in each actuator versus actuator number.

## *3.2. The Characteristics of a Composition of PSJA Array*

As shown in Figure 12, with the increase in the actuator number, the energy deposited in a single actuator decreases. As a result, the performance of the actuator ought to slow down. The history of jet velocity and gas density at the throat outlet is shown in Figure 13a,b, respectively. Obviously, the jet velocity and duration decreases with the increase in the actuator number. However, the gas density at the throat outlet increases with the increase of actuator number. As the gas kinetic energy is proportional to the gas density, the growth of gas density benefits the gas kinetic energy.

**Figure 13.** The history of gas velocity and gas density at the throat outlet versus actuator number. (**a**) is the gas velocity and (**b**) is the gas density.

The jet peak velocity, duration time, and mass ejected with different numbers of PSJA are calculated and results are as shown in Figure 14a. The peak velocity is defined as the maximum velocity at the throat outlet. The definition of the jet duration time can be found in Figure 8a. The mass ejected means the total mass of gas left in the PSJA chamber throughout the ejection process. These parameters decrease with the increase in the actuator number, but the rate at which they decline slows. These parameters are normalized by dividing the value with that of only one actuator, as shown in Figure 14b. This figure shows that the mass ejected is greatly influenced by the actuator. The influence of actuator number on jet duration time is weakest compared with the other two parameters. When the actuator increases to 20, the jet duration time decreases by 33%, while the peak velocity and mass ejected decrease by 54% and 62%, respectively.

**Figure 14.** The gas parameters versus actuator number. (**a**) is the real value and (**b**) is the normalized value.

The jet kinetic energy can be calculated as shown in Equation (25), where *ρ*o(*t*) and *u*(*t*) are the jet density and velocity at the throat outlet face; *A*<sup>t</sup> indicates the throat outlet area; *T*jet represents the jet duration time. The variation of jet kinetic energy of a single actuator with actuator numbers is shown in Figure 15. As it is proportional to the third power of the velocity, the kinetic energy decreases quickly with the increase in the actuator number. When the actuator number increases to 20, the jet kinetic energy decreases by 92%.

$$E\_{\rm kinetic} = \int\_{0}^{T\_{\rm jet}} 0.5 \rho\_{\rm o}(t) u^{3}(t) A\_{\rm t} dt \tag{25}$$

**Figure 15.** The gas kinetic energy versus actuator number.

## *3.3. The Energy Characteristics of PSJA Array*

The total kinetic energy of all actuators is calculated using Equation (26), which varies with the actuator number, as shown in Figure 16. In Equation (26), *k* is the concrete actuator; and *n* is the number of PSJA actuators in the PSJA array. Firstly, the total kinetic energy increases with the number of actuators. Then it begins to decrease, suggesting that the maximum total kinetic energy is not necessarily produced by the array containing more PSJAs and an optimal number of PSJAs exists. The characteristic of total kinetic energy is different from that of the total discharge energy, as shown in Figure 11, which will be investigated below.

**Figure 16.** The total gas kinetic energy versus actuator number.

In the whole working process of the actuator, the energy stored in the capacitor transforms to the kinetic energy indirectly, as shown in Figure 17. The thermodynamic recycle efficiency *η*<sup>t</sup> is defined in Equation (27) follows. The variation of *η*<sup>t</sup> with different actuator number is shown in Figure 18. As the actuator number increases, the energy deposited in an actuator decreases. Based on the research results of Zong et al. [26], the thermodynamic efficiency decreases with the decrease of heating energy. With the increase of actuator number, the heating energy decreases. Therefore, the thermodynamic efficiency decreases when there are more PSJAs. As a result, the total gas kinetic energy doesn't always increase with an increase of actuator number.

$$
\eta\_t = \frac{E\_{\text{total}}}{Q\_h} \mathcal{O} \tag{27}
$$

**Figure 17.** The total gas kinetic energy versus actuator number.

**Figure 18.** The thermodynamic efficiency versus actuator number.

#### *3.4. The Influence of Extra Resistance*

It is observed that the total extra resistance, including the wire resistance and the equivalent serial resistance of the discharge capacitor, affects the performance of the PSJA array greatly. In this part, these factors are studied. The extra resistance is set as 0.1 Ω and 10 Ω in case 1 and case 2 respectively.

The total discharge energy and discharge efficiency versus actuator number with different extra resistance are shown in Figure 19. In case 1, the extra resistance is only 0.1 Ω, and the discharge efficiency of one actuator is 72%. In contrast, the discharge efficiency of one actuator in case 2 (larger extra resistance) is as low as 5%. With the increase in the actuator number, the discharge efficiency in the two cases increases. However, the increased levels are different. When the actuator number increases to 20, the discharge efficiency in case 1 increases to 94%, which is 30% larger than the initial discharge efficiency (72%). But in case 2, the discharge efficiency increases to 38%, which is almost six times larger than the initial discharge efficiency (5%).

**Figure 19.** The total discharge energy and efficiency versus actuator number with different extra resistance.

The variation of total gas kinetic energy with actuator number under different extra resistance is shown in Figure 20. The revealed trends are completely different when different extra resistance is used. In case 1, the total gas kinetic energy decreases with the increase in the actuator number. However, in case 2, the opposite trend happens. Instead of decreasing, the total kinetic energy increases

when more actuators are used. In case 2, the total gas kinetic energy with 20 actuators is nearly four times of that with one actuator.

**Figure 20.** The total gas kinetic energy versus actuator number with different extra resistance.

The reason for the distinctive variations revealed in Figure 20 is given as follows. The discharge efficiency increases with actuator number, but the thermodynamic efficiency decreases, while the variation speed is different. The discharge efficiency and thermodynamics efficiency with different actuator numbers are normalized by that with only one actuator, as shown in Figure 21. The larger the extra resistance is, the larger the increase in the discharge efficiency is, and the smaller the decrease in the thermodynamic efficiency is. As a result, with a large extra resistance, the total efficiency grows with the increase in the actuator number. However, with a small extra resistance, the total efficiency decreases with the increase in the actuator number.

**Figure 21.** The normalized discharge efficiency and normalized thermodynamic efficiency versus actuator number with different extra resistance.

In conclusion, when the extra resistance is large, the energy efficiency of the PSJA array is higher than that of conventional PSJA. Actually, in the practical application, the extra resistance is relatively large. In the experiment on control of the shock boundary layer interaction using PSJA by Greene et al. [27], a resistance (max 1.5 kΩ) is adopted in the discharge circuit. In the experiment of Narayanaswamy et al. [26,28,29], a ballast resistance is also adopted in the discharge circuit. In the experiment of Jichul Shin [30,31], a ballast resistance (max several kilo Ohm) is adopted in the discharge circuit as well. Therefore, there is reason to believe the PSJA array driven by the multichannel discharge technology can play an important role in the practical flow control application.

## **4. Conclusions**

In this paper, based on the electronic theory, the conservation equations of mass, momentum, energy, and the lumped capacitance method, the multichannel discharge model, and plasma synthesis jet actuator model are put forward. Coupled with the two models, an analytical model to simulate the characteristics of the PSJA array is developed. Based on this model, the influence of actuator number in the PSJA array is investigated. The main conclusions are as follows.

With the increase in the actuator number, the discharge energy in a single actuator decreases. As a result, the jet strength induced by a single actuator decreases. When the actuator number increases from 1 to 20, the weakening extent of mass ejected, peak jet velocity, and jet duration time is 62%, 54%, and 33%, respectively. Therefore, the influence of the actuator number on the jet duration time is the weakest.

As the actuator increases, the discharge efficiency and thermodynamic efficiency show different tendencies. The discharge efficiency increases, while the thermodynamic efficiency decreases. As a result, the total energy efficiency doesn't always increase with the increase in the actuator number. However, in practical application, the extra resistance in the discharge circuit is large (to the order of kilo Ohms), and the discharge efficiency of one channel discharge is low. Therefore, it can be concluded that the PSJA array is more suitable to the practical application.

**Author Contributions:** Formal analysis, H.S.; Investigation, S.Z.; Supervision, Z.Z. and S.H.; Validation, W.Y. and L.Y.

**Funding:** This research was funded by the National Natural Science Foundation of China (91541120, 91641204, and 51790511) and Royal Society (IE150612).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Multibody Simulation for the Vibration Analysis of a Turbocharged Diesel Engine**

**Enrico Armentani 1,\*, Francesco Caputo 2, Luca Esposito 3, Venanzio Giannella <sup>4</sup> and Roberto Citarella <sup>4</sup>**


Received: 15 May 2018; Accepted: 17 July 2018; Published: 20 July 2018

**Abstract:** In this paper, a multibody calculation methodology has been applied to the vibration analysis of a 4-cylinder, 4-stroke, turbocharged diesel engine, with a simulation driven study of the angular speed variation of a crankshaft under consideration of different modeling assumptions. Moreover, time dependent simulation results, evaluated at the engine supports, are condensed to a vibration index and compared with experimental results, obtaining satisfactory outcomes. The modal analysis also considers the damping aspects and has been conducted using a multibody model created with the software AVL/EXCITE. The influence of crankshaft torsional frequencies on the rotational speed behavior has been evaluated in order to reduce the vibration phenomena. The focus of this work is related to industrial aspects since, for an existing and commercialized engine, a numerical and experimental complex study has been performed to enable design improvements aimed at reducing noise and vibrations. Existing procedures and algorithms are combined here to reach the abovementioned objectives in the most efficient way.

**Keywords:** vibration analysis; FEM; multibody simulations

## **1. Introduction**

Structural and acoustic modeling methods [1], used to predict the vibroacoustics and aeroacoustics performance of aerospace [2] and automotive systems, have become a key tool in the design process.

Car engine dynamics are explored especially in the low frequency range and at low engine speeds, where the direct vibration transmission by the engine mounts is a critical excitation mechanism, but raising the maximum analyzed frequency constitutes an important industrial challenge for automotive and aerospace industry.

There are several examples in the literature describing a combination of multi-body simulation, including flexible FE subsystems, to perform dynamics and acoustic simulation of mechanical components.

In [3], a multibody model of a valve train system for a car engine was developed in order to analyze its resonant vibrational behavior. The dynamics of the valve train system were analyzed, applying a ramp to the crankshaft with a variable engine speed and considering the motion unsmoothness deriving from the inertia and combustion pressure loads.

In [4] the development of a numerical model to predict the noise radiating from manual gearboxes due to gear rattle is provided. The measured data are used to identify and reproduce the input excitation that is primarily generated from engine combustion forces. The dynamic interaction of the gearbox components, including flywheel, input/output shafts, contacting gear-pairs, bearings, and flexible housing is modeled using flexible multibody techniques. The acoustic response to the vibration of the gearbox housing is then predicted using vibroacoustic techniques.

In [5,6] the numerical modeling of noise radiated by an engine and by a car body, respectively, using the so-called Acoustic Transfer Vectors and Modal Acoustic Transfer Vectors concepts is presented. The dynamics of the engine are described using a finite element model loaded with an RPM-dependent excitation.

In [7], a technique to compute the noise radiated by a large truck engine is presented: the vibration is computed using a commercial FEM code, whose results are subsequently used in an acoustic radiation computation.

In [8], a vibroacoustic numerical and experimental analysis was carried out for the chain cover of a low powered four-cylinder four-stroke diesel engine. A Boundary Element (BE) model of the chain cover was realized to determine the chain cover noise emission, starting from the previously calculated structural vibrations. The numerical vibroacoustic outcomes were compared with those observed experimentally, obtaining a good correlation.

In [9,10], a novel and effective mathematical model and advanced analytical approaches to achieve a more accurate prediction of the spiral bevel gear dynamic response were developed to investigate the underlying physics affecting gear mesh and gear dynamic response generation and transmissibility.

In [11], an in-depth investigation of the dynamical load sharing behaviors of a four-planetary gear system with multi-floating components was provided.

In [12], the impact of coupling engine structure with rotating components on the engine noise and vibrations across all rpm and frequency range is illustrated.

The need for accurate models forms a major obstacle to the implementation of cylinder balancing methods for engines with a high number of cylinders [13]. This is due to closely spaced cylinder firings and the fact that the crankshaft dynamics cannot be ignored, partly due to the increased length of the crankshaft, and partly because analysis of higher frequency components is required to obtain sufficient information for balancing the cylinder-wise torque contributions. This deformation can assume significant values depending on the engine-load configuration (load change, crankshaft stiffness, kind of aspiration of the engine), and as such it is of great importance for safe engine operation.

In [14], an experimentally validated diesel engine simulation code is used to study and evaluate the importance of a notable engine dynamic issue, i.e., the crankshaft torsional (angular) deformations during turbocharged diesel engine operation, owing to the difference between instantaneous engine and load (resistance) torques. The analysis aims ultimately in studying the phenomena under the very demanding, and often critical, transient operating conditions. Details are provided concerning the underlying mechanism of the crankshaft torsional deformations during steady state and transient operation.

In [15], a methodology for predicting the piston to liner contact in running engines by means of MBD (Multi-Body Dynamics) and FEM is presented. In addition to the mathematical modeling of the excitation, the paper describes the transfer mechanisms of the piston slap phenomenon. Thus, the model is extended in order to analyze vibration transfer via engine structure. Results of simulation work show structure surface velocity levels and their contribution to integral levels in different frequency bands.

In [16], the mathematical modeling of body structures and the calculation of the non-linear connecting forces resulting from elastohydrodynamic contacts between piston-liner and shaft-bearing is described. Results of parametric studies, e.g., the influence of piston surface profile on the contact mechanism between piston and liner, are shown. In [17], the coupling of elastic multibody simulations (including elastohydrodynamic interactions) with finite element based vibration and acoustic analysis is presented. In [18], the emphasis is on the integration of the kinetic reactions arising from the tribological conjunction of the dynamics of engine subsystems, piston, and crankshaft.

In [19], the CAE capabilities in the simulation of the dynamic and acoustic behavior of an engine, with a focus on the relative merits of modification and full-scale structural/acoustic optimization of engine, are presented.

In [20], an investigation of the diesel engine combustion-related fault detection capability of crankshaft torsional vibrations is presented: the torsional vibration amplitudes are used to superimpose the mass and gas torque; further mass and gas torque analysis is used to detect fault in the operating engine. The engine dynamics are analyzed with a focus on the low frequency range and low engine speeds, when vibration transmission through engine mounts becomes critical.

In [21,22], a detailed multi-body numerical model of an engine prototype is used to characterize the whole engine dynamic behavior in terms of forces and velocities. A combined usage of FEM and multi body methodologies is adopted for the dynamic analysis: both crankshaft and cylinder block are considered as flexible bodies, whereas all the other components are considered as rigid elements.

The focus of this work is related to industrial aspects because, for an existing and commercialized engine, a numerical and experimental complex study has been performed to enable design improvements aimed at reducing noise emission and vibrations. Existing procedures and algorithms are combined here to reach the abovementioned objectives in the most efficient way.

## **2. Problem Description and Modeling Approach**

The analyzed problem concerns a numerical study of the vibrations of an in-line 4-cylinder, 4-strokes, internal combustion turbocharged diesel engine, to be used as a first step for future vibroacoustic analyses. Instead of a complex numerical analysis of the system through a direct FEM, which could involve prohibitive number of degrees of freedom, the authors present a modeling approach exhibiting a gradually increasing complexity. Starting from rigid body analysis and introducing progressively the elastic behavior of the various subsystems, using then a modal FEM analysis, the impact on accuracy of such modeling refinements comes out. The originality of the overall approach is related to the combination of various numerical approaches with an experimental analysis, for a very complex system.

The analyses are focused on the 2nd, 4th, 6th, and 8th order of motion irregularities that are analyzed at the flywheel and pulley. In addition, the time-dependent simulation results at the engine support brackets (engine bracket, gearbox bracket and differential bracket) are evaluated, condensed to a vibration index and finally compared with the experimental results.

A Multi-Body Dynamic Simulation (MBDS) of the crank train was used to characterize its dynamic behavior, starting from engine geometrical data and the available combustion loads, with both mechanical and combustion forces acting simultaneously on the crankshaft.

In order to examine the vibration behavior of the considered internal combustion engine, the behavior of the crankshaft was specifically analyzed with a focus on its torsional vibrations, namely calculating the motion irregularities of the crankshaft itself.

For the issue at hand, the modeled components were: crankshaft, pistons, connecting rods, main and connecting rod bearings, engine mounts (also named "support brackets"), contact stiffness between piston and cylinder.

The first realized model assumed rigid bodies, allowing assessing the course of the motion irregularities of the crankshaft at low frequencies. With a maximum regime for the engine established at 4500 rpm, we considered 0–300 Hz a low frequency range: component resonances in this range can be activated by the 2nd and 4th engine order. The considered mid-frequency range was 300–700 Hz: resonances in this range can mainly be activated by the 6th and 8th engine order. High frequencies, higher than 700 Hz, can only be activated by aeroacoustics phenomena.

The next step was the introduction of a crankshaft flexible model, leveraging on a FEM approach, in order to evaluate more accurately its dynamic behavior at higher frequencies. The frequency range of

interest was anyway limited to the low-medium frequencies because the whole engine was modeled as rigid with the only exception of the crankshaft. Focusing on the low-medium frequency range, it was reasonable to assume a negligible influence from the modal behavior of other components than the crankshaft.

Subsequently, the crankshaft FE model was further refined taking into consideration the presence of another component, the clutch: this allowed obtaining more accurate support brackets vibrations.

Finally, a numerical-experimental correlation for the validation of the numerical model was carried out.

Hypermesh and Abaqus [23] codes were respectively used for the FE modeling and modal analysis, whereas the realization of the multibody (MB) model and the forced vibrations analysis was demanded to AVL/Excite code [24].

The MB code considered the behavior of individual bodies as linear elastic; such bodies can be subject to both large rigid body motions and small deformations.

The applied external forces come from the pressure cycle of the combustion gases [25]: from Figure 1 it is possible to appreciate the cycle variations at different regimes. All the forces of an inertial nature are calculated internally by the code according to the actual speeds and accelerations of the bodies. Calculations performed for each operating regime and in the time domain provided the displacement, speed, and acceleration time histories of all the points of the system.

**Figure 1.** Experimental pressure cycles at different regimes for the engine under analysis.

#### *2.1. System Dynamic Reduction*

The adopted flexible MB approach is 'floating-frame-of-reference component-mode-synthesis (FFR-CMS)', whose theoretical description can be found in the MB books by Shabana [26,27]. The Craig Bampton dynamic condensation [28] was adopted. The first 50 eigenforms are used for the modal reduction (such number is judged sufficient for the scope).

The free-free natural frequencies of the reduced model are calculated and provided as input to the MB analysis: when the flexible body is connected, through the various joints, to the other bodies included in the overall model, the constrained natural frequencies are calculated and used to solve the forced analysis by a modal response approach.

The concept underlying the creation of a MB model is to subdivide a mechanical system, having an overall nonlinear elastic behavior, into linear elastic sub-systems and to concentrate nonlinearities in the connections between them. The elastic bodies are represented by the condensed matrices of the corresponding FE model. A number of nodes, usually those chosen in condensation, having a mass and connected to each other by massless springs and dampers, discretizes each elastic body.

In Figure 2a,b the schematization of the constraints and components making up the engine are shown, with highlight of engine mounts ("FTAB" body) and corresponding location with respect to the engine (modeled as a rigid body "RI3D"). The condensation nodes of the flexible crankshaft ("CON6") are shown in Figure 2c.

**Figure 2.** *Cont*.

**Figure 2.** Graphic representation of bodies and connections (**a**,**b**) and highlight of crankshaft points used as master nodes for the reduction process (**c**).

For the damping matrix, the MB software refers to the so-called 'proportional damping', known as "Rayleigh damping" (Figure 3). This can be defined by using two frequencies and the damping ratio at these frequencies. The function of the damping ratio ε is as follows:

$$\mathfrak{E}\_{\mathsf{i}} = (\mathsf{a} + \mathsf{b} \ \omega\_{\mathsf{i}} \mathsf{2}) / (2 \ \omega\_{\mathsf{i}}) \text{ with } \omega\_{\mathsf{i}} = 2 \ \mathsf{\pi} \mathsf{f}\_{\mathsf{i}}.$$

where


**Figure 3.** Rayleigh damping vs. frequency (Hz) for the crankshaft.

## *2.2. Models of Increasing Complexity*

For the analysis of the vibration behavior of the internal combustion engine under examination, three models of increasing complexity were realized:


**Figure 4.** Crank train with highlight of force introduction.

**Figure 5.** *Cont*.

**Figure 5.** (**a**) Powertrain with highlight of support positions: (**b**) engine bracket; (**c**) gearbox bracket; (**d**) differential bracket.

**Figure 6.** Highlight of added point N.2, geometrically coincident with the N.1 but rigidly connected to the seismic mass.

For all the models, the remaining powertrain components were considered as rigid, with assigned characteristics such as weight, mass center position and inertia. The gearbox was not modeled, neither was present in the experimental bench test (an engine brake was adopted).

Therefore, the purpose of model SOL1 is purely methodological, while models SOL2 and SOL3 are those on which the analyses of interest and outcomes comparison were carried out.

## 2.2.1. Model SOL1

For model SOL1, some components were schematized by an appropriate rigid body, such as a crankshaft, pulley, flywheel, clutch, connecting rods, pistons, piston pins and rings, or power unit (base, head, oil pan, distribution components, and accessories). In these cases, the specifications of geometric characteristics, values of inertia, mass, and position of the centers of gravity were sufficient.

Moreover, the chassis was schematized by means of a "Generic Body" that represents the car chassis and can be considered as a constraint on the ground to which the powertrain must be clamped by three brackets: the engine, gearbox, and differential brackets (Figure 2a,b). In particular, "Generic Body" is a tool of EXCITE customized for the specific application by introducing the proper parameter settings.

In order to represent the constraints between adjacent bodies, appropriate joints were used (Figure 2a,b). In the case of bearings, they were modeled by a revolution joint, with a given spring stiffness (related to the bearing stiffness) and damping value (related to the oil film characteristics). The stiffness values of the bearings were available from AVL database and selected for this specific engine version based on the following default values:


The joint force in radial direction is determined in the plane perpendicular to the main rotational axis as a function of the previously described parameters.

The crankshaft was connected to the crankcase by means of "Revolute Joints" that simulate the presence of the main bearings; this type of constraint allows us to recreate a linear contact between the two bodies connected by a spring-damper modeling. In particular, the value of the clearance between shaft and crankcase and the stiffness and damping of the spring-damper system were inserted: these values depend on the maximum pressure in the combustion chamber, the bore, and the maximum force discharged on the single main bearings.

No modeling of the oil film was provided, namely, hydrodynamics of the fluid in the main bearings (e.g., between piston and cylinder liner and the connecting rod) was not considered. Actually, in the literature it is proposed to take into account even the elastohydrodynamic interactions for some applications, because, for vibroacoustic analysis, some local modes, like those in correspondence of the bearings, can play an important role especially when considering the crankshaft bending behavior. Anyway, the focus of this work was solely on the shaft torsional behavior. In fact, the elastohydrodynamic interactions generally do not influence the torsional behavior because, due to a sufficiently high crankshaft stiffness, there is no coupling between the bending modes affected by elastohydrodynamic interactions and the torsional modes. Moreover, elastohydrodynamic allowance would cause a huge increment (three times) of run times.

The bearings used in the engines to counteract the axial thrusts generated during motion were schematized with the constraint "Axial Thrust Bearing"; also for this case, clearance, stiffness, and damping values were provided. However, axial thrusts are considerably lower than those deriving from the vertical thrust, generated in the combustion chamber and acting on the main bearings (with this constraint, it was possible to connect a node of the shaft to more nodes of the crankcase).

The connecting rods were connected to the power unit through two constraints: the first component was a prismatic guide schematized through the Piston-Liner Guidance constraint and used to simulate the constraint between rods and plungers; the second component simulated the connecting rod bearing, for which the joint used for main bearings was taken into consideration.

The presented modeling does not calculate the torsional critical speeds, since the motion of this calculation model is rigid, and therefore it is a zero-pulse model. Therefore, the elements making up the system do not deform and rotate at the same speed.

## 2.2.2. Model SOL2

Model SOL2 introduces the flexibility of the crankshaft, getting a different dynamic behavior of the system and allowing the evaluation of those critical torsional speeds non-visible from model SOL1. This element of flexibility was introduced into the analysis through an FE modeling of the body crankshaft. The crankshaft comprised 150,000 nodes, 6 degrees of freedom (DOFs) per node, whereas the reduced condensed model comprised 50 nodes and 300 DOFs.

The free-free natural frequencies of the reduced model were calculated and provided as input for the MB analysis. When the flexible body was connected through the various joints to the other bodies included in the model, the constrained natural frequencies were calculated and used to solve the modal frequency response analysis.

In such a model (Figure 7), the clutch unit was modeled as a concentrated mass. Figure 8 shows the FE models of the flywheel and pulley.

The adopted mesh utilizes tetrahedral elements with an average size equal to 0.5 mm.

**Figure 7.** *Cont*.

**Figure 7.** FE model of the shaft (**a**) with clutch mechanism modeled as a concentrated mass (**b**).

**Figure 8.** FE model of the flywheel (**a**) and pulley (a cutout is shown) (**b**).

## 2.2.3. Model SOL3

Model SOL3 comprises the FE modeling of the clutch mechanism (Figure 9) in addition to the crankshaft already present in the SOL2 model.

Crankshaft and clutch mechanism meshes comprised nearly 450,000 tetrahedral and hexahedral elements and nearly 140,000 nodes. The salient features of the three models, SOL1, SOL2 and SOL3, are also condensed in Table 1.

**Figure 9.** FE model of the shaft with discretized clutch mechanism.


#### **Table 1.** Characteristics of the calculation models.

### *2.3. Experimental Setup*

The experimental measurements were performed on a test bench (Figure 10), considering in particular a measurement point located on the powertrain flywheel. Therefore, in order to have a reliable comparison, the experimental motion irregularities were compared with the numerical ones in a reference node located in the mass center of the flywheel (node 108 shown in Figure 6). Calibration between numerical and experimental tests was based on the measured torque.

A preliminary check on the correlation of free-free numerical and experimental crankshaft eigenmodes was performed with satisfactory outcomes.

**Figure 10.** Test bench for the engine under analysis.

## **3. Results**

## *3.1. Model SOL1*

Although not leading to useful outcomes regarding the determination of the critical torsional speeds, Model SOL1 allowed us to obtain important outcomes for the model verification purposes, i.e., motion irregularity of crankshaft and global vibration index for the brackets.

It is well known that in an internal combustion engine, the engine torque generated by gas and inertia forces is periodic and subject to breakdown into a Fourier series as a sum of harmonics that can be more or less relevant to the engine dynamics according to the associated amplitude, phase, and order value.

## 3.1.1. Motion Irregularities

The outcomes of the simulation obtained with a rigid crankshaft (SOL1), focus on the study of the motion irregularities detected in the center of gravity of flywheel (node 108), hub (node 1), and seismic mass (node 2), as in Figure 6.

The even orders are the most relevant for the engine torque of a 4-stroke, 4-cylinder engine; therefore, their influence should be considered for motion irregularities assessment. Consequently, an evaluation of the magnitude of the 2nd engine order of crankshaft motion irregularity vs. speed was primarily carried out for the engine under test (Figure 11a).

The Campbell diagram represents the frequency spectrum of a non-stationary signal, in which the frequency is on the abscissa and the average angular velocity is on the ordinate, whereas the plotted data indicate the oscillation amplitude of the variable under test (in this case Δrpm). Campbell diagram for the motion irregularity evaluated on the flywheel (Figure 11b) clearly shows the prevalence of the 2nd engine order. It is interesting to observe that the 6th order seems to play a negligible role in the motion irregularities; however, this will be invalidated when the allowance for modal behavior is included in the analysis, triggering resonance with a torsional mode.

**Figure 11.** Motion irregularities—model SOL1: (**a**) 2nd engine order with highlight of Δrpm vs. speed (rpm) and (**b**) Campbell diagram of the motion irregularities with highlight of 2nd, 4th and 6th engine orders.

## 3.1.2. Global Vibration Index

A global vibration index, named as weighted sum, was also calculated for the brackets and defined as follows:

$$X(f) = \sqrt{\sum\_{i=1}^{n} w\_i^2(f) \cdot x\_i^2(f)},\tag{1}$$

where:

*wi*(*f*) is the weighting factor for a given direction and position,

*xi*(*f*) is the rms acceleration for a given direction and position,

*n* is the total number of positions and directions (*x*, *y* and *z*) = number of brackets × 3.

The weighting factors used in Equation (1) depend on the vehicle in which the powertrain has to be installed and are generally provided by the vehicle platform unit. Without such information, as in the present case, the coefficients are all set to 1 (uniform weighing of the structural transmission paths).

The second order of the global vibration index of each bracket (engine, gearbox and differential brackets) is shown in Figure 12, with a contribution that is nearly frequency independent for gearbox and differential brackets, but increasing with frequency when considering the engine bracket. The corresponding global vibration index is shown in Figure 13 with reference to the amplitudes of numerical and experimental acceleration irregularities, and therefore to the amount of the transmitted vibration. Such a good numerical–experimental agreement can be explained by observing that the 2◦ engine order cannot trigger a torsional mode and consequently, when considering the second order effect, there is no loss of accuracy by neglecting the crankshaft flexibility.

**Figure 12.** *Cont*.

**Figure 12.** Numerical bracket vibration indexes *X*(*f*) [m/s2] for the second engine order: (**a**) gear box; (**b**) engine; (**c**) differential.

**Figure 13.** Numerical–experimental comparison of the brackets global vibration index for the second engine order.

Model SOL1 allowed us to verify some operating parameters, such as the motion irregularities and the global vibration index for the brackets; it also allowed us to highlight that, as expected, the second engine order harmonic force is predominant when neglecting the flexibility of any component.

## *3.2. Model SOL2*

## 3.2.1. Motion Irregularities

From the modal analysis, it is possible to observe that, for the engine under testing in its operating frequency range, considering the real constraints (as enabled by the multibody analysis), the first two torsional modes are at 305 and 545 Hz, respectively (Figure 14). In the former, the hub and the seismic mass oscillate in phase, whereas their oscillation is out of phase for the latter. Anyway, the differences between free-free and constrained torsional eigenfrequencies are very low (nearly 15 Hz considering the first torsional mode with vertical and transversal springs at the bearings' locations); on the contrary, a relevant difference exists between the free-free and constrained bending modes.

**Figure 14.** 1st (**a**) and 2nd (**b**) torsional eigenmodes at 305 and 545 Hz, respectively.

The primary objective here is to detect how the amplitude of the motion irregularities is affected by the shaft torsional eigenfrequencies: in fact, in correspondence of resonances, the irregularities may have values so high as to make them critical to the vibroacoustic behavior of the whole powertrain.

The outcomes of model SOL2 are compared with the test data to assess the accuracy of this modeling. Considering that the experimental measurements are only available on the flywheel center, the numerical–experimental comparison can only involve node 108 (Figure 6).

The 2nd order of motion irregularity does not display any peak in the range 1250–5000 rpm (Figure 15); this is expected since the maximum involved excitation frequency,

$$f = 2 \times 5000/60 = 166.7 \text{ Hz}, \pi$$

is not sufficiently high to trigger the first torsional eigenmode at 305 Hz; consequently, SOL2 provided the same results of SOL1 and a satisfactory correlation with experimental measurements made at flywheel (node 108), as in Figure 16. The minimum of the motion irregularities is at nearly 4250 rpm: in correspondence with such a regime, the combustion forces and inertia forces are almost in equilibrium.

**Figure 15.** Numerical motion irregularities for the 2nd engine order at nodes 1, 2, and 108, as provided by model SOL2.

**Figure 16.** Numerical-experimental comparison of the motion irregularities for the 2nd engine order at node 108.

On the contrary, the 4th order of motion irregularity does display a peak in the range 1250–5000 rpm at nearly 4500 rpm (Figure 16), corresponding with the frequency

$$f = 4 \times 4500/60 = 300 \text{ Hz}, \pi$$

which is sufficiently close to the 1st torsional eigenfrequency (305 Hz) to confirm the hypothesis of a resonant behavior for the crankshaft. Such a peak only manifests with reference to nodes 1 and 2, but is negligible for node 108 (Figure 17).

A good correspondence between the numerical and experimental results is provided by SOL2 (Figure 18) for regimes higher than 2000 rpm, even if a slight underestimation of the irregularities is obtained. The lack of numerical vs. experimental correlation at low regimes is expected due to the turbolag phenomenon. In fact, during the bench test based on an acceleration from the min to the max regime, the turbo did not properly work at low regime (lower than 2000 rpm), whereas the simulation cannot allow for such behavior since it is based on a quasi-stationary variation of loading conditions.

It is worth noting that, as expected, the amplitude of oscillation for node 2 is higher than that of node 1: in fact, hub and seismic mass oscillate in phase, with the latter being more peripheral and consequently having larger oscillations.

**Figure 17.** Numerical motion irregularities for the 4th engine order at nodes 1, 2, and 108, as provided by model SOL2.

**Figure 18.** Numerical-experimental comparison of the motion irregularities for the 4th engine order at node 108, as provided by model SOL2.

Exploring the 6th order, a numerical peak at nearly 3050 rpm can be observed [18], corresponding with the frequency

$$f = 6 \times \text{3050/60} = \text{305 Hz}, \text{m}$$

which, again, is coincident with the 1st torsional eigenfrequency (305 Hz), confirming the hypothesis of resonant behavior for the crankshaft.

The experimental peak is at 3100 rpm, whereas the numerical peak is at 3050 rpm (Figure 19): this discrepancy can be attributed to the "discrete" acquisition procedure adopted in the test.

The discrepancies of motion irregularities at the low regimes can be surely explained by the turbolag phenomena, whereas the non-negligible underestimation provided by the simulation at high regimes (Figure 20) might be caused by the fact that the only crankshaft is modeled as flexible.

**Figure 19.** Numerical motion irregularities for the 6th engine order at nodes 1, 2, and 108, as provided by model SOL2.

**Figure 20.** Numerical-experimental comparison of the motion irregularities for the 6th engine order at node 108, as provided by model SOL2.

Regarding the 8th order, two peaks at nearly 2250 and 4250 rpm can be found (Figure 22), corresponding with the frequencies

$$f\_1 = 8 \times 2250 / 60 = 300 \text{ Hz},$$

$$f\_2 = 8 \times 4250 / 60 = 567 \text{ Hz}.$$

These frequencies are sufficiently close to the 1st (305 Hz) and 2nd (545 Hz) torsional eigenfrequencies, confirming the hypothesis of resonant behavior for the crankshaft.

Again, a slight underestimation of motion irregularities provided by the simulation is evident (Figure 21).

In conclusion, the 2nd torsional eigenfrequency cannot be triggered by orders lower than the 8th; moreover, orders higher than the 8th are neglected because their amplitude is sufficiently low to prevent a relevant impact when triggering the torsional eigenfrequencies.

**Figure 21.** Numerical-experimental comparison of the motion irregularities for the 8th engine order at node 108, as provided by model SOL2.

**Figure 22.** Numerical motion irregularities for the 8th engine order at nodes 1, 2, and 108, as provided by model SOL2.

## *3.3. Model SOL3*

Regarding model SOL3, the outcomes of the 2nd engine orders motion irregularities do not significantly differ from those observed for model SOL2. In contrast, different results between the two models are obtained for the 4th engine order (Figure 23). In particular, a different trend of the brackets vibrations due to the explicit FE modeling of the clutch mechanism is observed.

**Figure 23.** Global vibration index for the 4th engine order brackets: comparison between models SOL2 and SOL3.

## **4. Conclusions**

In this work, three internal combustion turbocharged diesel engine multibody models of increasing complexity were developed and their vibration behavior observed.

The development of these multibody models allowed us to analyze the support brackets vibrations and crankshaft motion irregularities in order to proceed to a design optimization solution in production and/or the development phase.

The choice of adding the only flexibility of crankshaft, flywheel, and pulley is related to the range of frequencies of interest (up to 650 Hz) and to the exclusive focus on the impact of the crankshaft torsional behavior (the flexibility of other powertrain components comes out at higher frequencies). Consequently, in the frequency range of interest, it was not necessary to adopt complex models for the calculation of both the vibrational index at the engine brackets and the irregularity motion of crankshaft, with a substantial reduction in the computational burden.

The obtained outcomes of the three analyzed models can be summarized as follows:


The next steps for this work would be allowance for the flexibility of engine block and engine subsystems in order to calculate the impact of its resonances on the monitored vibrational parameters at higher frequencies. Clearly, such improvements would also require further considerations of the damping factors to be used in the calculations, to be fine-tuned by numerical experimental correlations performed on the single components.

**Author Contributions:** Conceptualization: E.A. and R.C.; Data curation: E.A., F.C., L.E., V.G. and R.C.; Formal Analysis: E.A. and R.C.; Investigation: E.A. and R.C.; Methodology: E.A. and R.C.; Resources: E.A.; Project administration: E.A.; Validation: E.A., F.C., L.E., V.G. and R.C.; Visualization: E.A., F.C., L.E. and V.G.; Writing—original draft: E.A.; Funding acquisition: L.E.; Writing—review and editing: V.G. and R.C.; Supervision: R.C.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Study on Optimal Compensation Design for Meteorological Satellites in the Presence of Periodic Disturbance**

## **Shijie Xu 1, Naigang Cui 2, Youhua Fan 1,\* and Yingzi Guan <sup>2</sup>**


Received: 5 June 2018; Accepted: 18 July 2018; Published: 20 July 2018

## **Featured Application: This work may be used to improve the accuracy of attitude stability of flexible meteorological satellites.**

**Abstract:** Periodic disturbance may cause serious effects on the attitude of meteorological satellites, and the attenuation of periodic disturbance is required. In this paper, a fundamental study on the optimal design of constant compensations against known-law periodic disturbance for meteorological satellites is investigated. An analytical solution for the relationship between the frequency and amplitude ratios and the response of a typical second-order vibration system is firstly derived. The compensation and disturbance torques are determined according to practical engineering. The criteria for designing the optimal compensations are based on the analytical and simulation results. Then, the criteria are applied to a flexible spacecraft actuated by constant control torque in the presence of sustained periodic disturbance. The optimal compensation torque parameters for spacecraft are acquired based on former criteria. The compensate effectiveness of the optimal compensation torque is provided and compared with results of other selections in the frequency and amplitude ratio domain. Numerical simulation results and experimental results clearly demonstrate the good performance of proposed criteria. This work provides a significant reference for the vibration attenuation of meteorological satellites in the present of periodic disturbance.

**Keywords:** flexible spacecraft; periodic disturbance compensation; compensate torque design; vibration attenuation; reaction wheel.

## **1. Introduction**

Meteorological satellites play a vital role in our daily life. To send back high-quality images, satellites have to meet stringent demands related to attitude accuracy [1–3]. However, periodic disturbances from spaceborne equipment result in difficulty in retaining stability in the attitude of spacecraft. The application of spaceborne microwave imaging instruments has substantially improved the image definition, but poor effects on attitude accuracy are also brought about by their known-law rotating imbalance [4–6]. Furthermore, characteristics like flexibility and low weight make spacecraft vulnerable to vibrations caused by numerous sources [7]. In this paper, we focus on the compensation of known-law long-lasting periodic disturbances and attenuation of vibrations to improve the attitude accuracy of meteorological satellites.

Many studies have been done on the attitude control of meteorological satellites. Mitigation of the periodic disturbances by compensating torques has been proposed by researchers over the course of the past years. Zhang et al. [8] presented a fault tolerant control for external disturbances and actuator failures. This controller is the combination and application of the free weighting matrices method and LMI (linear matrix inequality) technique. However, this simply designed controller has not been well verified for flexible spacecraft. Cong et al. [9] presented an improved sliding mode attitude control algorithm with an exponential time-varying sliding surface. This controller enhanced the accuracy of maneuver in the presence of disturbances and parametric uncertainties. To reduce the dependence of estimate accuracy on a disturbance model, Hu et al. [10] designed a sliding mode disturbance observer for the attitude control of a rigid spacecraft in the presence of disturbance and actuator saturation. The adaptive controller has also been adopted. Unfortunately, the studies neglected the flexibility of the structure. The attitude control and disturbance rejection have been studied for flexible spacecraft during attitude maneuvering by Zhong et al. [11]. An improved compensator is developed to achieve asymptotic disturbance rejection. Then, a new robust state feedback controller is presented for the attitude maneuver. The proposed strategy achieved the maneuver process in the presence of small amplitude disturbances. Zhang et al. [12] proposed a generalized-disturbance rejection (GDR) control for vibration suppression of flexible structures in the presence of unknown continuous disturbances. The unknown disturbances are estimated by a generalized-disturbance rejection which is defined by a refined state space model. Chen et al. [13] focused on the disturbance observer-based control method to attenuate or reject disturbances for a spacecraft attitude control system. The disturbance observer is used to estimate the disturbances. Also, the feedforward compensation is given according to the corresponding results. Ma et al. [14] developed an adaptive failure compensation approach for the attitude control of satellites in the presence of disturbance and uncertain failures of actuators. This research mainly focuses on the disturbances from actuator failures; external disturbances have not been sufficiently considered. Lau et al. [15] improved the research on disturbance identification and rejection. The dipole-type disturbance rejection filter is adopted and applied in [15]. Subsequently, two simple close-loop system identification methods are introduced to improve the robustness of the rejection filter relative to frequency uncertainty. Chen et al. [16] made contributions to the rejection of general multi-periodic disturbances. The multi-periodic disturbances are estimated by a new presented method. Based on the estimation results, a new adaptive control method is applied for disturbance rejection. The residual vibration can be controlled in a specified range under the provided strategy. However, this approach assumes that all periods of disturbances should be known, which is difficult in practical engineering. Although the preceding strategies have made some achievements in disturbance rejection, parameter magnitudes considered in research, like the variation range of disturbance frequency and disturbance amplitudes, are not large enough. Furthermore, corresponding experiments verifying their effectiveness are absent in most literature, and it is unpractical to output complicated desired compensating torques by any actuators.

Maneuver actuators may also affect the attitude accuracy. Speaking of the actuators, jet thrusters are inappropriate for persistent disturbance. Reaction wheels are extensively used due to their advantage in generating continuous and precise torques without expending the propellant [17,18]. Although it is easy for torque mode reaction wheels to output a constant value torque, evident faults exist when tracking an expected time varying torque [19]. Most of the present papers failed to consider the capability of reaction wheels. These complex control output torques may be unfulfillable in common practical engineering. On the other hand, problems like time delay may be solved by the using of flywheel actuators. Sabatini et al. [20,21] tackled the problem of controlling time delayed systems. A prediction of state at the time of the control application has been performed to compensate for the time delay. A Kalman filter and the GNC loop were used to realize this purpose. Sabatini et al. also conducted research to test and improve the robustness of presented algorithm. Simulations and experiments were provided to confirm the effectiveness. Those works may cover the shortage caused by the time delay. Further, for meteorological satellites in the presence of known laws of periodic disturbances, those methods are inappropriate and superfluous. Based on those facts, for the compensation of periodic disturbed torques, constant reaction wheel torques may be the most practical and appropriate way. Unfortunately, work is limited on this issue in the literature.

Another way to improve the attitude accuracy of meteorological satellites is to actively suppress the residual vibrations of flexible appendages. This kind of strategy has been widely used in the attitude maneuver of spacecraft. Numerous studies have been presented for vibration suppression of flexible appendages. Smart materials, such as Magnetorheological Material (MR) [22], Piezoelectric Material (PM) [23,24], have been investigated and fabricated to eliminate the residual vibrations during the past decades. A wide range of control schemes has been proposed for the use of smart materials, such as positive position feedback (PPF) control [25], sliding mode control (SMC) [26], H-infinity control [27], and so on. However, these actuators and sensors have changed the structural property by mechanical interference, which is undesirable in many practical engineering applications. An active vibration suppression method called component synthesis vibration suppression (CSVS) can suppress vibration without additional structures attached to flexible structures. This method was first proposed by Liu et al. [28]. Then, Hu et al. [29] extended the application of the CSVS method in the attitude maneuver of spacecraft. However, as a feedforward control scheme, the application of the CSVS method is limited to the sensitivity to parameter uncertainties. Similar to the CSVS method, the input shaping technique is also used for vibration suppression. Kong et al. [30] extended the input shaping method by combining it with feedback control for the vibration control of flexible spacecraft during attitude maneuver. Their work remedied the shortcoming of input shaping as it is a feedforward control scheme. However, this technique is suitable for zero to zero maneuver. It is not appropriate for problems forced in this paper in which the disturbance is nonstop. All those tactics mentioned above have application in vibration attenuation of flexible appendages, but their application conditions are restrictive. They are not suitable for situations in which satellites are actuated by sustained disturbances. Motivated by the preceding discussion, finding an appropriate continuous compensating torque to attenuate the vibration of flexible appendages is the optimal option. However, direct studies about relationships between selection of compensate torques and vibration attenuation effectiveness for meteorological satellites are scarce.

In this paper, based on the requirement of meteorological satellites and practical limitations, step torques given by reaction wheels are adopted to compensate for the sinusoidal periodic disturbance. Two coefficients, frequency ratio and amplitude ratio, are introduced as parameters. A typical 2-order vibration system and a flexible spacecraft system are introduced. The optimal compensate torques and corresponding effect will be discussed by analytic methods. Residual vibrations are regarded as the standard of the disturbance compensate effect. To verify the effectiveness of the proposed strategy, simulation results are presented. We aim to determine the relationship between the selection of the two ratios and disturbance compensate effectiveness. Furthermore, an experimental study is provided to verify the effectiveness of the proposed strategy. It is hoped that the work will be helpful for solving practical engineering problems. In addition, preliminary investigations in this paper have limitations, and further discussions will be summarized in our next study.

The remaining part of this paper is organized as follows. Section 2 gives the mathematical modelling. Analytical relationships between ratios and residual vibration are given in Section 3. Numerical simulations, experiment results and discussions are presented in Section 3. Section 4 concludes the paper.

## **2. Materials and Methods**

#### *2.1. Modelling of the Flexible Spacrcraft*

Three-axis spacecrafts with flexible appendages have a complex model. It is difficult for us to obtain the relationship between vibrations and control torques. To stress the key point of this paper, here we provided a simplified spacecraft model: A flywheel actuated single-axis gas suspending rotary table (SGSRT) is applied. The diagram of the SGSRT is shown in Figure 1.

**Figure 1.** Diagram of SGSRT.

As is shown in Figure 1, this system is mainly consists of a rigid hub and a flexible beam. A rigid beam is fixed on the opposite side of flexible beam. An eccentric mass is rotating around the tip of the rigid beam at a constant angular velocity. A reaction wheel is placed at the centre of the hub. Define *OXYZ* as the inertial frame and *oxy* as the frame fixed on the edge of rigid hub. *θ*(*t*) is the rotation angle of the hub. *y*(*x*, *t*) is the displacement of place *x* in frame *oxy*. Ω and Ω*<sup>m</sup>* are the angular velocity of rigid hub and eccentric mass, respectively. The flexible beam is assumed to be an Euler-Bernoulli beam. The Newton-Euler method is used to build the mathematical model of the system. The dynamic model of the flexible system shown in Figure 1 is

$$\begin{cases} \begin{aligned} I\dot{\boldsymbol{\Omega}} + \int \boldsymbol{\dot{r}\_b} \times \boldsymbol{a}\_b d\boldsymbol{x} + \boldsymbol{r}\_{\rm m} \times \boldsymbol{a}\_{\rm m} \boldsymbol{m}\_m &= \boldsymbol{T} \\\ E I \frac{\partial^4 \boldsymbol{y}}{\partial \boldsymbol{x}^4} + \sigma \left[ \frac{d^2 \boldsymbol{y}}{dt^2} + (R\_{hub} + \boldsymbol{x}) \frac{d^2 \boldsymbol{\theta}}{dt^2} \right] &= 0 \end{aligned} \end{cases} \tag{1}$$

where *J* is the rotational inertia matrix of the rigid hub. *rb* is the vector from the centre of the hub to an arbitrary mass point in the flexible beam. *ab* is the acceleration of mass points. *r<sup>m</sup>* is the vector from the centre of the hub to the eccentric mass. *a<sup>m</sup>* and *mm* are the acceleration and mass of the eccentric mass, respectively. *T* is the actuator torque. *Rhub* is the radius of the rigid hub. *EI* is the flexural rigidity of the beam and *σ* is the surface density

As the SGSRT is adopted in this study, Equation (1) can be simplified because it is a single-axis rotation system. The simplified model can be expressed as:

$$\begin{cases} \begin{array}{c} J\ddot{\theta} + \int (R\_{hub} + \mathbf{x}) \ddot{y} dm + \mathcal{W}\_m \Omega\_m^2 = T\\ \, \, \, EI \frac{\partial^4 y}{\partial x^4} + \sigma \left[ \frac{d^2 y}{dt^2} + (R\_{hub} + \mathbf{x}) \frac{d^2 \theta}{dt^2} \right] = 0 \end{array} \tag{2}$$

where *J* is the rotation inertia of the system relative to the axis of rotation. *Wm* represents the influence on the system caused by eccentric mass. It can be expressed as *Wm* = *Pmσ<sup>m</sup>* sin(*ωmt*)*mm*. *Pm* is a constant parameter about the position of the eccentric mass. *σ<sup>m</sup>* is the horizontal distance from the eccentric mass to the tip of the rigid beam. *y* = *y*(*x*, *t*) is the displacement of the flexible beam.

Using the assumed-mode method, the displacement of the flexible beam can be expressed as:

$$\mathbf{y} = [D]\boldsymbol{\eta} \tag{3}$$

where [*D*] is a matrix consisting of modes of the flexible beam, *η* is the modal coordinate vector. In consideration of the orthogonality of the mode. Equation (2) can be shown as:

$$\begin{cases} \begin{array}{c} J\ddot{\boldsymbol{\theta}} + \mathcal{S}\ddot{\boldsymbol{\eta}} + \mathcal{W}\_{m}\boldsymbol{\Omega}\_{m}^{2} = \boldsymbol{T} \\ \boldsymbol{\delta}^{T}\ddot{\boldsymbol{\theta}} + \ddot{\boldsymbol{\eta}} + \mathcal{C}\dot{\boldsymbol{\eta}} + \mathbf{K}\boldsymbol{\eta} = \boldsymbol{0} \end{array} \tag{4}$$

where *δ* is the coupling matrix of rotation and vibration. *C* and *K* are damping and stiffness matrices, respectively. It needs to be emphasized that the reaction wheel is assumed to be ideal equipment. That means it can generate ideal square torque. Furthermore, the influence on the system caused by eccentric mass is in sinusoidal law according to its expression. Hence, the two known-law kinds of torque will be adopted in the residual part.

#### *2.2. Modelling of the Typical Second Order Vibration System*

A typical second order vibration system is also used due to its advantage in distinct mathematical analytical solution. It can be shown as:

$$m\ddot{\mathbf{x}} + c\dot{\mathbf{x}} + k\mathbf{x} = F \tag{5}$$

where *m*, *c* and *k* are the respective mass, damping and stiffness of the system. *F* is the mechanical load.

#### **3. Results and Discussion**

## *3.1. Analysis and Design of Optimal Compensate Criterion*

For many studies about vibration suppression of flexible spacecraft, their simplified mathematical model can be expressed as a 2-order vibration system, for example, like the model used in [7]. This means the 2-order vibration system can be used to illustrate the vibration characters of flexible spacecraft. In this section, the typical second order vibration system will be used for analysis of the effectiveness of disturbance compensation and to confirm the optimal compensate criterion. Sinusoidal disturbance and square compensation are adopted based on practical engineering. It should be noted that the frequencies of two kinds of torques are assumed to be the same.

Assume there are two periodic exciting forces, *F*<sup>1</sup> and *F*2, acting on the system. Then the system shown in Equation (5) can be expressed as:

$$
\ddot{x}\,m\ddot{x} + c\dot{x} + kx = F\_1 + F\_2\tag{6}
$$

where

$$F\_1 = \begin{cases} B, & \frac{2j\pi}{\omega\_f} \le t < \frac{(2j+1)\pi}{\omega\_f} \\ -B, & \frac{(2j+1)\pi}{\omega\_f} \le t < \frac{2(j+1)\pi}{\omega\_f} \end{cases} \text{ } B \ge 0, \ j = 0, 1, 2, \dots$$

$$F\_2 = -A \sin\left(\omega\_f t\right), \text{ } A \ge 0$$

*ω<sup>f</sup>* is the natural frequency of both *F*<sup>1</sup> and *F*2. *B* and *A* are the amplitudes of the square force and sinusoidal force, respectively. Periodic square force *F*<sup>1</sup> can be expressed by Fourier series as:

$$F\_1 = \frac{a\_0}{2} + \sum\_{n=1}^{\infty} \left( a\_n \cos n\omega\_f t + b\_n \sin n\omega\_f t \right) \tag{7}$$

where the coefficients are:

$$a\_0 = 0, \ a\_n = 0, \ b\_n = \begin{cases} \frac{4B}{n\pi}, & n = 1, 3, 5, \dots \\ 0 & n = 2, 4, 6, \dots \end{cases}$$

In this series, the fundamental harmonic has the maximal amplitude. Then, the right side of the equal sign in Equation (6) can be expressed as:

$$F\_1 + F\_2 = \left(\frac{4B}{\pi} - A\right)\sin\omega\_f t + \frac{4B}{\pi}\left(\frac{1}{3}\sin 3\omega\_f t + \frac{1}{5}\sin 5\omega\_f t + \dotsb\right) \tag{8}$$

Here, we introduce an auxiliary parameter called amplitude ratio, which can be expressed as *β* = *<sup>B</sup> <sup>A</sup>* . Substitute Equation (8) into Equation (6), and replace *B* by *β*, we have:

$$x\ddot{\mathbf{x}} + c\dot{\mathbf{x}} + k\mathbf{x} = \sum\_{i=1}^{\infty} P\_i \sin\left((2i-1)\omega\_f t\right) \tag{9}$$

where

$$P\_i = \begin{cases} \frac{(4\beta - \pi)A}{\pi}, & i = 1\\ \frac{4A\beta}{(2i - 1)\pi}, & i = 2, 3, 4, \dots \end{cases}$$

Solving Equation (9), we obtain

$$\mathbf{x}(t) = \sum\_{i=1}^{\infty} \left\{ \begin{aligned} \varepsilon^{-\frac{\tau}{\zeta\_i}\omega\_{\mathrm{tr}}t} & \left( \mathbf{x}\_0 \cos \omega\_d t + \frac{\dot{\mathbf{x}}\_0 + \frac{\tau}{\zeta\_d} \omega\_{\mathrm{tr}} \mathbf{x}\_0}{\omega\_d} \sin \omega\_d t \right) \\ & + Q\_i \varepsilon^{-\frac{\tau}{\zeta\_i}\omega\_{\mathrm{tr}}t} \left[ \sin \varrho\_i \cos \omega\_d t + \frac{\omega\_{\mathrm{tr}}}{\omega\_d} (\zeta\_i \sin \varrho\_i - \lambda (2i - 1) \cos \varrho\_i) \sin \omega\_d t \right] \\ & + Q\_i \sin \left( (2i - 1)\omega\_f t - \varrho\_i \right) \end{aligned} \tag{10}$$

where *x*<sup>0</sup> = *x*(0), . *<sup>x</sup>*<sup>0</sup> <sup>=</sup> . *<sup>x</sup>*(0) are the initial conditions. *<sup>λ</sup>* <sup>=</sup> *<sup>ω</sup><sup>f</sup> <sup>ω</sup><sup>n</sup>* is the frequency ratio.

$$Q\_i = \frac{P\_i/k}{\sqrt{\left(1 - \left(2i - 1\right)^2 \lambda^2\right)^2 + \left((4i - 2)\zeta\lambda\right)^2}}, \ q\_i = \arctan\frac{(4i - 2)\zeta\lambda}{1 - \left(2i - 1\right)^2 \lambda^2}$$

For zero initial conditions system, the solution shown in Equation (10) can be expressed as:

$$\mathbf{x}(t) = \sum\_{i=1}^{\infty} \left\{ \begin{array}{l} Q\_i e^{-\int\_{\omega} \omega\_i t} \left[ \sin \varrho\_i \cos \omega\_d t + \frac{\omega\_0}{\omega\_d} (\zeta \sin \varrho\_i - \lambda (2i - 1) \cos \varrho\_i) \sin \omega\_d t \right] \\ + Q\_i \sin \left( (2i - 1) \omega\_f t - \varrho\_i \right) \end{array} \right\} \tag{11}$$

Equation (11) is the final displacement formula. The first line of Equation (11) will be decreased because of damping, and only the second line will remain when system becomes stable. We can find the displacement is strongly associated with the frequency and amplitude ratio, and the fundamental harmonic is vital in the compensated force design. Numerical simulations were provided to provide the laws.

'Simulation' in 'Matlab' was used here for exact residual vibration under different combinations of the frequency ratio and amplitude ratio. For the system of Equation (6), we set its natural frequency and damping ratio as *<sup>ω</sup>sp* <sup>=</sup> <sup>√</sup>0.9 and *<sup>ζ</sup>sp* <sup>=</sup> 0.01/*ωsp*, where the subscript 'sp' means mass-spring system. The amplitude value of sinusoidal disturbance *A* is set to be 1. In this way, only amplitude ratio and frequency ratio are required. It should be noted that the parameter used here is idealized. The detailed relevance of criterion design and damping will not be considered too much here. We may study those problems in the future. The compensation effect may not notable or excessive if the two ratios are too small or too big. The range of the amplitude ratio is selected from 0.3 to 1 while the range of the frequency ratio is selected from 0.3 to 3. The tip vibration is chosen to be the standard to show the effectiveness of disturbance torque suppression. Simulation results are exhibited in Figures 2–5.

Figure 2a,b are the total displacement results in two representations. The variation tendency is clear. At frequency ratio *λ* = 0.3 and *λ* = 1, residual vibrations sharply increased. These are resonances caused by different excitations. The system is mainly receiving two kinds of excitations. One is the composition of the fundamental harmonic of square force and another is the disturbance. Both of the two excitations have

the same frequency *ω<sup>f</sup>* . Resonance appeared when *λ* = 1. Another kind of excitation is the remaining harmonics of square force. Frequency of the second order harmonic almost meets the fundamental natural frequency of system when *λ* = 0.3. As a result, resonance will appear again in the frequency ratio-domain. For each frequency ratio, there is always an amplitude ratio where the residual vibration reaches the minimum. For most frequency ratios, the best compensation effect appeared near the amplitude ratio *β* = 0.7. The result for frequency ratio *λ* = 1 is a typical case. For every frequency ratio, the fundamental harmonic of the square force compensates the sinusoidal force near *β* = 0.7. Those facts illustrated that *β* = 0.7 is the optimal option for most feasible frequency ratios.

**Figure 2.** Displacement amplitude results in two forms: (**a**) In surface representation; (**b**) In marker representation.

Figure 3a,b provides the accelerations with different combinations of the amplitude and frequency ratio. From the results, we can find that the variation laws are the same as those for the displacement results. Resonances are obvious at *λ* = 1 and *λ* = 0.3. Beyond that, values vary homogeneously at different frequency ratios. The condition under which residual vibrations always have a minimum value for each frequency ratio is also distinct, especially for frequency ratios larger than 1. The minimum values are at about the centre of the amplitude ratio range. To show the results more clearly, Figures 4 and 5 show the comparisons between different ratios.

**Figure 3.** Acceleration amplitude results in two forms: (**a**) In surface representation; (**b**) In marker representation.

Figure 4a,b are parts of the results shown for the frequency ratio-domain. Four equally distributed amplitude ratios are selected. Data in two figures have the same variation laws. Situations when amplitude ratio *β* = 0.7 are the best as their vibrations are the smallest for almost all the frequency ratios. That is because the sinusoidal disturbance is properly counteracted by the fundamental harmonic of the step force. Resonance is conspicuous at frequency ratio *λ* = 0.3 and *λ* = 1.

**Figure 4.** Results for the frequency ratio-domain: (**a**) Displacement amplitude results; (**b**) Acceleration amplitude results.

**Figure 5.** Results for the amplitude ratio-domain: (**a**) Displacement amplitude results; (**b**) Acceleration amplitude results.

Figure 5a,b are parts of the results shown for the amplitude ratio-domain. Each curve has its own minimum, but not all the minimums are located near *β* = 0.7. For frequency ratios greater than 1, the minimum appeared near *β* = 0.7. For frequency ratios less than 1, the minimum amplitude ratios are not fixed. Floating windows in Figure 5a,b illustrate that for some frequency ratios, their minimum may appear near *β* = 0.2. However, the corresponding frequency ratios do not frequently occur in practical engineering.

Results indicate that there is always an optimal combination of the amplitude and frequency ratio for disturbance compensation. In general, we can list the optimal compensate criterion: The optimal amplitude ratio is located within 0.6 to 0.7 for most of the frequency ratios. Besides, the frequency ratio should avoid 0.3 and 1, where it may cause resonance.

## *3.2. Simulations for Flexible Spacecraft*

The former subsection demonstrated that the square compensation force can be used to mitigate the vibration caused by sinusoidal disturbance. To demonstrate the effectiveness of compensation in flexible spacecrafts, numerical simulations and experiments are implemented. In this subsection, the flexible spacecraft model in Equation (4) is adopted, and simulation results are presented.

The system is actuated by a sinusoidal disturbance and a square compensate torque produced by a periodic rotated eccentric mass and a reaction wheel, respectively. Residual vibrations with different combinations of frequency and amplitude ratios will be given to show the effectiveness.

Here, the torque *T* can be expressed as:

$$T = T\_1 + T\_2 \tag{12}$$

where *T*<sup>1</sup> and *T*<sup>2</sup> have the same expression as *F*<sup>1</sup> and *F*<sup>2</sup> in Equation (4), respectively. *T*<sup>1</sup> is given by the ideal reaction wheel. *T*<sup>2</sup> is from the periodic eccentric mass. *T*<sup>1</sup> and *T*<sup>2</sup> can be expressed as:

$$T\_1 = \begin{cases} \mathsf{U}, & \frac{2j\pi}{\omega\_m} \le t < \frac{(2j+1)\pi}{\omega\_m} \\ -\mathsf{U}, & \frac{(2j+1)\pi}{\omega\_m} \le t < \frac{2(j+1)\pi}{\omega\_m} \\ & T\_2 = -\mathsf{W}\_m \Omega\_m^2 \end{cases} \\ \mathsf{U} \ge 0, \ j = 0, 1, 2, \dots$$

Equation (4) can be further expressed as:

$$\begin{cases} \begin{array}{c} J\ddot{\boldsymbol{\theta}} + \boldsymbol{\delta\ddot{\boldsymbol{\eta}}} = \sum\_{i=1}^{\infty} P\_i \sin((2i-1)\omega\_m \boldsymbol{t}) \\ \boldsymbol{\delta\boldsymbol{T}} \ddot{\boldsymbol{\theta}} + \ddot{\boldsymbol{\eta}} + \mathbf{C} \dot{\boldsymbol{\eta}} + \mathbf{K} \boldsymbol{\eta} = \mathbf{0} \end{array} \end{cases} \tag{13}$$

where

$$P\_i = \begin{cases} \frac{(4\beta - \pi)}{\frac{\pi}{T}} \Omega\_m^2 P\_m \sigma\_m m\_{m\nu} & i = 1\\ \frac{4\Omega\_m^2 P\_m \sigma\_m m\_m \beta}{(2i - 1)\pi}, & i = 2, 3, 4, \dots\\ \beta = \frac{U}{\Omega\_m^2 P\_m \sigma\_m m\_m} \end{cases}$$

Substitute the first part of Equation (13) into the second formula of Equation (13), we will have:

$$\left(I\_3 - \frac{\delta^T \delta}{l}\right)\ddot{\eta} + \mathbf{C}\dot{\eta} + \mathbf{K}\eta = -\frac{\delta^T}{l}\sum\_{i=1}^{\infty} P\_i \sin((2i-1)\omega\_m t) \tag{14}$$

Modal solutions are complex to be predicted analytically, not to mention the multi-dimension modal. Thus, we will directly present the simulation results. Parts of parameters are measured values of experimental equipment used in our study. Others are calculated in a mathematical way. The moment of inertia of the rigid hub is *<sup>J</sup>* = 15 kg·m2. The damping ratios of the vibration are *<sup>ζ</sup>*<sup>1</sup> = *<sup>ζ</sup>*<sup>2</sup> = *<sup>ζ</sup>*<sup>3</sup> = 0.05, the natural frequencies are *ωn*<sup>1</sup> = 2.11 rad/s, *ωn*<sup>2</sup> = 13.28 rad/s, *ωn*<sup>3</sup> = 37.15 rad/s. The coupling matrix is *δ* = 2.65 0.54 0.21 kg·m2. In simulations, the mass of rotating eccentric mass is *mm* = 0.4 kg, and its eccentric radius is *σ<sup>m</sup>* = 0.08 m. The states of system are initialised to be zero.

The maximum of sinusoidal disturbance will vary in a certain range depending on each frequency ratio. Square compensation torques were assumed to be ideally outputted, which means it can perfectly fit the sinusoidal disturbance. The frequency ratio range is set from 0.3 to 3.2, and the amplitude ratio range is set from 0.3 to 1. Those settings are based on the results in the previous subsection. Amplitudes of flexible beam tip residual vibration under different combinations of frequency and amplitude ratios are calculated by 'Simulation' in 'Matlab' software. Simulation results are illustrated in Figures 6 and 7.

Figures 6a,b show the displacement of the flexible beam tip with all the combinations of frequency and amplitude ratio in this simulation. Similar to the results displayed in the former subsection, when the amplitude ratio is near *β* = 0.7, the displacement of the flexible tip is minimal. This law is particularly obvious when the frequency ratio exceeds 0.5. Those results just proved the optimal criterion. Variation amplitudes are so weak that the trend is not so clear when the frequency ratio is smaller than 0.5.

**Figure 6.** Displacement amplitude results in two forms: (**a**) In surface representation; (**b**) In marker representation.

Figure 7a,b are the amplitude history of the flexible beam tip acceleration. The tip acceleration is a real effect when the amplitude ratio is located near *β* = 0.7. The frequency of the system is fixed, so the increasing frequency ratio heightened the rotation velocity of the eccentric mass, and also the amplitude of torques. As a result, values of accelerations increase with increasing frequency ratio.

**Figure 7.** Acceleration amplitude results in two forms: (**a**) In surface representation; (**b**) In marker representation.

Figure 8a,b show the amplitudes of residual vibration in the frequency ratio-domain with selected amplitude ratios.

**Figure 8.** Results in the frequency ratio-domain: (**a**) Displacement amplitude results; (**b**) Acceleration amplitude results.

The variation tendencies in Figure 8a,b are analogous to those in Figure 4a,b. Situations with *β* = 0.7 in the frequency ratio-domain are much smaller than those for the other three selected amplitude ratios. Resonance is obvious at *λ* = 1 for all the amplitude ratios even though the amplitude of the synthetic torque is small. The fundamental harmonic of square torque is adopted to compensate the sinusoidal disturbance while the other orders of the harmonic become the new disturbances. Frequency of the second order harmonic meets the fundamental natural frequency of the system precisely, but amplitudes of remaining orders of the harmonic are so small to cause evident vibrations. As a result, resonances at *λ* = 0.3, which are shown in floating windows, are not obvious. Although higher Fourier series result in new disturbances, the main effect is from the first order.

Figure 9a,b shows us how amplitudes vary in the amplitude ratio-domain under special frequency ratios.

**Figure 9.** Results for the amplitude ratio-domain: (**a**) Displacement amplitude results; (**b**) Acceleration amplitude results.

Results revealed by Figure 9a,b demonstrated the relationship between the amplitude ratio and the level of residual vibration. The results all have the feature that the displacement amplitude always has a minimum in the amplitude ratio-domain. According to data in Figure 9a, *β* = 0.7 is an optimal region for the compensation of disturbances. The optimal amplitude ratio for tip acceleration, shown in Figure 9b, appears to be somewhat different, but the best combination of the amplitude and frequency ratio does exist despite that the values are different for frequency ratios. All the simulation results demonstrated this for the most frequently occurring situations.

### *3.3. Experimental Results*

Experiments were performed on an SGSRT to verify the proposed criterion and the simulation results. The equipment is shown in Figure 10. This systems mainly consists of a single-axis rotation hub, a flexible aluminium beam (2480 × 200 × 2.5 mm), a torque mode reaction wheel (*T*max = 0.065 Nm) whose rotation axis is coincident with the hub, a periodic rotation dish with an eccentric mass to create the disturbance, and equipment on the hub for measuring, control and data exchange.

**Figure 10.** Experimental setup of a single-axis rotary table.

The hub is suspended by a gas bearing. A zero-gravity environment is required, and the table is limited to rotate on the horizontal plane. The residual vibration is measured by an accelerometer fixed on the tip of a flexible beam. The model of the accelerometer (±2 *g*,0 − 400 *Hz*) is 2220-020, produced by Silicon Designe Inc. (Kirkland, WA, USA). Eccentric mass is chosen to be 0.4 kg and 0.25 kg for different frequency ratios, and the periodic rotation dish is placed 0.6 m from the rotation axis of the system. This is the best setting to generate sinusoidal disturbance whose amplitude is within the limits of the reaction wheel.

Frequencies of the experimental equipment have a close relationship with the experimental results. To identify the frequencies of the experimental equipment, free vibration tests were performed, and the fast Fourier transform (FFT) was applied. Figure 11a,b shows the measured tip acceleration data and the FFT result.

**Figure 11.** Tip free vibration and FFT analysis: (**a**) Measured tip free vibration; (**b**) FFT of tip free vibration.

Results show that the first two vibration orders are obvious. The measured and the calculated frequencies are shown in Table 1.


**Table 1.** The first and second frequencies (rad/s).

The errors of the first two frequencies between measured and calculated results are 2.781% and 2.19%, respectively. In view of the difference between the experimental equipment and mathematical model, these errors are acceptable and negligible.

According to the previous results, there are 2 typical varying laws in the amplitude-ratio domain. For situations *λ* < 0.5, variation is not obvious. This is because the torques are too small for flexible systems. Similarly, the torques are too large when *λ* is too high. Considering those situations and the output ability of the experimental equipment, three frequency ratios 1.3, 0.9 and 0.6 were selected for the experiment. Thus, three cases of the experiment corresponding to three selected frequency ratios were provided. On the other hand, simulations illustrated that the optimal amplitude ratio for most frequency ratios is near 0.7. Thus, three amplitude ratios near 0.7 were also selected for experimental cases.

Simulations were performed with real parameters in the experiment. Zero initial conditions were provided for all the experiments. Square compensation torques from the reaction wheel and sinusoidal disturbance were started at the same time. The measured accelerations of the flexible beam tip are shown in Figures 12–14.

## 3.3.1. Case A

In this case, *λ* = 0.6, *σ<sup>m</sup>* = 0.08 m, *mm* = 0.4 kg. Three experiments were performed with *β* = 0.6, 0.8, 1, respectively. The tip residual accelerations are shown in Figure 12.

In case A, the frequency ratio is selected as 0.6. In this case, three typical experiments are conducted with three amplitude ratios. The measured tip accelerations are illustrated in Figure 12. The tip accelerations were recorded. It can be seen from Figure 12 that the vibration amplitude when *β* = 0.6 is the smallest. Vibration amplitude when *β* = 0.8 is slightly greater than that when *β* = 0.6. However, when the amplitude ratio is 1, the residual vibrations are much more obvious than before.

**Figure 12.** Measured accelerations of the flexible beam tip in Case A.

## 3.3.2. Case B

In this case, *λ* = 0.9, *σ<sup>m</sup>* = 0.04 m, *mm* = 0.4 kg. Three experiments were performed with *β* = 0.6, 0.8, 1, respectively. The tip residual accelerations are shown in Figure 13.

**Figure 13.** Measured accelerations of the flexible beam tip in Case B.

In case B, the envelope curves are obvious. That is because the frequency of torques is close to the natural frequency of the system. As the amplitude ratio increases, the maxima of the residual vibration increase.

## 3.3.3. Case C

In this case, *λ* = 1.3, *σ<sup>m</sup>* = 0.04 m, *mm* = 0.25 kg. Three experiments were performed with *β* = 0.6, 0.8, 1, respectively. The tip residual accelerations are shown in Figure 14.

In this case, the situations when *β* = 0.6 and *β* = 0.8 are almost the same, but when *β* = 1, the residual vibration amplitude increased evidently.

Experimental results illustrate the square torques are efficient to compensate sinusoidal disturbance, but the effects vary with different conditions. According to the results of experiments, for all the three cases, the compensation effects are almost the same when *β* = 0.6 and *β* = 0.8. Compared with that, situations when *β* = 1 are worse. The amplitudes of residual vibration during compensation are listed in Table 2.

**Figure 14.** Measured accelerations of the flexible beam tip in Case C.


**Table 2.** Tip acceleration amplitude with different frequency and amplitude ratios m/s2.

Recall the added special simulation results which corresponded to the experiments. Figure 15 presents the results in the amplitude ratio-domain. The point marked 'ex' means the experimental data. The experimental results were largely in line with the simulation results. Simulation results show that the maximum amplitude of residual vibration always has a minimum point for each frequency ratio. Despite that the number of experiments is limited, the effectiveness and veracity of the proposed criterion are proved. Both experiments and simulations have verified the veracity of our study. However, it should be noted that all the experiment results are from the SGSRT, and the limitation of single axis rotation can be avoided. This limitation will be solved in future studies.

**Figure 15.** Experimental results and corresponding simulation results.

## **4. Conclusions**

In this article, the optimal criterion of constant compensations has been investigated for flexible spacecraft in the presence of known-law sustained periodic disturbance and square torques. This criterion is mainly based on the mathematical analytical results of the typical vibration system. First, the constant square compensation torque is obtained using the Fourier expansion method, and the relationship between the response of the system and input torques is obtained. The optimal criterion is obtained according to the analysis and simulation results. The criterion exhibits advantages in a flexible spacecraft system in the presence of sustained periodic disturbance. In general, the optimal amplitude ratio is near 0.7 for most frequencies. The frequency ratio should avoid the points near 0.3 and 1. Both numerical simulations and experiments have been performed to verify the veracity of the obtained criterion. However, defects such as the disturbance being too simple and the SGSRT is equipment with only one axis rotation are obvious in our study, and extension studies will be investigated for those limitations in our future research.

**Author Contributions:** S.X. performed the simulations, analyzed the data and wrote the paper; N.C. directed the mathematical analysis; Y.F. directed the simulation analysis and supervised this study. Y.G. directed the experiments.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number [91648201].

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18