**Research on Leakage Location of Spacecraft in Orbit Based on Frequency Weighting Matrix Beamforming Algorithm by Lamb Waves**

**Lei Qi 1,2, Zhoumo Zeng 1, Yu Zhang 1,\*, Lichen Sun 2, Xiaobo Rui 1, Xin Li 3, Lina Wang 2, Tao Liu <sup>2</sup> and Guixuan Yue <sup>1</sup>**


Received: 13 January 2020; Accepted: 7 February 2020; Published: 11 February 2020

## **Featured Application: The Frequency Weighting Matrix Beamforming Algorithm proposed in this paper has the potential to locate leakages in large structures such as spacecraft, and the algorithm can adapt the characteristics of leaks.**

**Abstract:** Clashes between space debris and spacecraft in orbit may cause air leakages, which pose a substantial danger to the crew and the spacecraft. Lamb wave dispersion in spacecraft structures and the randomness of leak holes are the difficulties in leak location. To solve these problems, a frequency weighting matrix beamforming algorithm is proposed in this paper. The elastic Lamb waves that are caused by leakages are acquired by an 'L' shaped sensor array consisting of eight acoustic emission sensors. The angle of a leak can be obtained through the superposition of different time delays, and the intersection of two angles can be used to find the location of the leak. Traditional beamforming is improved by matching the wave speeds in different frequency bands and weightings according to the energy distribution. Narrowband filtering is used to delay overlay different signal speeds with different frequency bands via a dispersion curve. The weighting method is used to compensate the frequency band response of different leak holes. The detailed location algorithm process is introduced and verified by experiments. For 1.5 and 2 mm leak holes, location direction accuracies of 1.33◦ and 1.93◦ for one sensor array were obtained, respectively.

**Keywords:** leakage location; Lamb wave; beamforming; spacecraft in orbit

## **1. Introduction**

Orbital objects, such as meteoroids and space debris, are among the serious and inevitable threats to spacecraft [1]. Technically, collisions caused by debris larger than 10 cm can be avoided by using databases that are compiled by debris-tracking systems [2]. Meanwhile, strikes by objects that are smaller than 10 cm lack an efficient detecting approach. A collision may damage the bulkhead of the spacecraft, causing internal pressure to leak, which seriously threatens the flight and the safety of astronauts and leads to severe consequences [3]. Thus, the detection and location of leakages would be of significant use for spacecraft in orbit.

For the problem of leak source location, some related findings have been reported. According to different principles, the currently representative methods include the optical imaging method [4], the optical fiber method [5], the helium mass spectrometry method [6], and the acoustic detection method [7]. Compared with other methods, the acoustic detection method can achieve the non-destructive detection and real-time location of leakages in a larger range with a high location accuracy and stability. Moreover, acoustic sensors are easy to integrate with spacecraft structures, thus eliminating the need for complex processes such as laying optical fibers or using noble gases. Therefore, the acoustic detection method is a potential leakage location technology for spacecraft in orbit [8].

For leakage location based on acoustic methods, the key difficulty lies in the continuity of the leakage acoustic signal. For a collision with a time point of occurrence, location can be performed by calculating the time difference of arrival. To solve the above problems, Holland et al. [9,10] carried out research on cross-correlation and two-dimensional Fourier transform techniques to determine air leakages on the International Space Station (ISS). In the study, two 16 × 16 sensor arrays were applied to locate a leak. All possible cross-correlations between the 256 sensor positions in each array were recorded. A large number of sensors and millions of data samples led to a time- and finance-consuming task. Al-Jumaili et al. [11] proposed a T-mapping method. This method obtains time-domain eigenvalues from a prior meshing experiment and uses a clustering algorithm to identify the position of the signal source. This scheme has good adaptability to complex structures, but the initial workload is huge and it has no application potential for huge spacecraft.

In order to reduce the number of sensors and the algorithm complexity of the leak location scheme, the beamforming technique was firstly applied to the method of leakage detections with acoustic emission sensors by McLaskey et al. [12]. The beamforming algorithm has been primarily used in radar, sonar and exploratory seismology [13]. However, the specimen tested in their experiment was a steel-reinforced concrete bridge ramp instead of a thin plate. The work was based on Rayleigh wave theory, which is not applicable to the thin plates used in spacecraft. Tian et al. [14] tested a thin plate that was equipped with near-field beamforming analysis, which was suitable for Lamb wave theory. Burst acoustic emission signals that were created by breaking a mechanical pencil lead on the surface of the specimens were applied in their tests, and the propagation speeds of the signals were obtained by time difference of arrival (TDOA) techniques. This method is not ideal for a leak-generated continuous signal. However, the difficulties of using beamforming techniques for acoustic emission signal processing lie in the confirmation of wave speeds due to their dispersion. Zhang et al. [15] considered the effect of dispersion, but their scheme did not consider the diversity of the frequency bands that are generated by leakage in practical applications. Thus, it can be seen that the difficulties in the leakage location of thin plates such as those used spacecraft are the effects of dispersion on wave speed and the band adaptability of practical applications.

In order to solve the above problems and to improve the accuracy of the beamforming algorithm for leak location, a frequency weighting matrix beamforming (FWMB) algorithm is proposed in this paper. In this method, narrowband filtering is used to match the dispersion curve to achieve the matrix superposition of wave speeds at different frequencies. An energy evaluation on the frequency band of the signal is performed to determine the selected frequency band and its weighting coefficient matrix, and the weighting method is used to compensate for the frequency band response that is caused by different leaks to achieve a frequency band adaptation. The FWMB algorithm was validated by laboratory tests that used an L-shaped acoustic emission array, and its performance was evaluated. The method in this paper has the potential to be used in spacecraft in orbit due to its light weight and easy integration.

The remainder of this paper is organized as follows. In Section 2, the FWMB process is proposed. In Section 3, the experiment platform in the laboratory is introduced. In Section 4, the conducted experiments are described in order to demonstrate the performance of the FWMB algorithm. The conclusions are presented in Section 5.

## **2. Method**

## *2.1. Delay-and-Sum Beamforming Algorithm*

Beamforming is a mature digital signal processing technology that is commonly used in the fields of radar, sonar, and seismic detection. The basic idea of the beamforming algorithm is to calculate the energy of a wave in a certain hypothetical direction by summing the delays of the signal of each array element to give an estimate of the wave arrival direction. The beamforming method uses a sensor array with a fixed spatial position to measure the spatial sound field, and it processes the signals that are measured by each sensor to obtain detailed sound source information. In the leak detection method that is covered in this paper, the direction of a leak can be obtained by a set of sensors. The final location point can be obtained by the intersection of two sets of sensor arrays.

In beamforming algorithms, the choice of sensor array shape is critical. According to the characteristics of spacecraft structures, a planar array can be used for detection. The shape of a planar array is usually circular, square, triangular, cross, and the like. For spacecraft, the load needs to be as small and light as possible. Research by Cui et al. [16] showed that L-type arrays can obtain the best experimental results with the smallest number of sensors. Therefore, in the following research, an L-shaped array with 8 sensors was used as an example. Moreover, the method proposed in this paper is applicable to other array formats. The demonstration of the beamforming algorithm is shown in Figure 1.

**Figure 1.** Demonstration of beamforming.

As shown in Figure 1, the sensor array in the figure is labeled n, including #0–7. Among them, the #0 sensor is for reference. When a leak occurs, sound waves propagate through the plate and are acquired by the sensor. *f*(*t*,*n*) is defined as the signal that is acquired by each sensor at the time of *t*. The relative angle θ between the leak source and the array is defined as the actual direction of the leak source, and θ*'* represents the hypothetical direction.

When the distance between the sensor array and the leakage source is far greater than the size of the array, the signal waves are considered to propagate along parallel paths. Under this assumption, the signal sampled by the reference sensor #0 is the replica of time-delayed signals that are recorded by other array elements. The time delays between array elements are determined by both by the relative

distance from the signal source to array element #*n* compared to the reference sensor #0 under the direction of arrival and by the wave speeds. Usually, the direction of arrival is unknown; a hypothetical direction of arrival θ*'* is set. According to the geometric relationship of Figure 1, the time delay of each signal for #0 can be expressed as:

$$
\Delta t\_n = \frac{d\_{\rm tr}}{v} \,\prime\tag{1}
$$

where *dn* represents the relative distance as shown in Figure 1 and *v* represents the wave speed. *dn* can be defined as:

$$d\_n = \left\{ \begin{array}{c} n \times d \times \cos \theta', [n = 1, 2, 3] \\ (n - 3) \times d \times \sin \theta', [n = 4, 5, 6, 7] \end{array} \right. \tag{2}$$

where d represents the distance between the sensors; this distance is the same as the diameter of the sensor and is 8 mm.

If θ*'* is consistent with the true angle θ, the signal will be concentrated through the superposition of time delay. Therefore, the signals of each sensor are delayed and superimposed to obtain the following delayed sum signal *g*(*t*,θ*'*):

$$\log(t, \theta') = \sum\_{n=1}^{7} f(t - \Delta t, n) + f(t, 0), \tag{3}$$

The energy of the signal can be obtained by squaring and integrating the superimposed signal in the time domain. By scanning calculations for different hypothetical angles, an energy function *B* (θ*'*) related to the angle can be obtained:

$$B(\theta') = \int g^2(t, \theta') \, dt,\tag{4}$$

After obtaining the angle of the peak, it can be used as the result. The location of the leak source can be obtained by the intersection of the two arrays' straight lines.

#### *2.2. Dispersion of Lamb Wave*

Considering that the bulkhead structure of an orbiting spacecraft is a thin metal plate, due to its large curvature, it is approximated as a metal flat plate in this paper to simplify the research. In thin metal plates, ultrasonic waves propagate in the form of Lamb waves. Lamb waves have less attenuation than body waves, which is beneficial to the location of leaks. However, the propagation velocity of Lamb wave varies with wave frequency, which is called the dispersion characteristic of Lamb waves. Frequency dispersion complicates the characteristics of ultrasonic waves, and it is difficult to directly determine the speed of waves to complete the beamforming algorithm. As a consequence, studying the dispersion characteristics of Lamb waves is critical. There are two modes of waves—symmetric and anti-symmetric—that independently propagate in plates. By equipping a numerical method, equations known as the Rayleigh–Lamb frequency relations for both symmetric and anti-symmetric waves expressed as Equations (5) and (6) can be solved [17].

$$\text{Symmetric mode}: \frac{\tan(qh)}{q} + \frac{4k^2 p \tan(ph)}{\left(q^2 - k^2\right)^2} = 0,\tag{5}$$

$$\text{Anti-symmetric mode}:\ q\tan(qh) + \frac{\left(q^2 - k^2\right)^2 \tan(ph)}{4k^2p} = 0,\tag{6}$$

where *k* represents the wave number, ω is the angular frequency, and *d* is the thickness of the plate.

When the thickness of a plate and the frequency of the signals are given, a set of numerical solutions, denoted by *ka* and *ks*, are acquired by solving the Rayleigh–Lamb equations. *ka* and *ks* represent the wave numbers of anti-symmetric and symmetric waves, respectively. This phenomenon proves that multiple Lamb wave modes exist in the plate, and each mode has its own phase velocities and group velocities. The number of the numerical solutions increases with the frequency of the signals and the thickness of the plates, which means that the number of wave modes increases. Usually, symmetric modes are represented by S0, S1, ... , Sn, and anti-symmetric modes are represented by A0, A1, ... , An. In each mode, the corresponding phase velocity also varies with the frequency, and the phenomenon is called the frequency dispersion. Dispersion curves consist of a set of curves can be drawn to depict these characteristics. From the dispersion curves, the theoretical values of phase speeds at any particular frequency–thickness product of any specific mode can be obtained.

#### *2.3. Frequency Weighting Matrix Beamforming Algorithm*

The traditional beamforming algorithm and the dispersion phenomenon in the plate are introduced above. From the analysis, it can be understood that the traditional beamforming method has the following two difficulties for the leakage location of thin plates such as spacecraft shells:


The FWMB algorithm that is proposed in this paper solves the above problems through the following two angles:


A flowchart comparison between the traditional beamforming [18] and the FWMB algorithm is shown in Figure 2. In the figure, the main differences between the two algorithms are marked by a red border. The following describes the FWMB algorithm based on the traditional beamforming algorithm with some improvements in detail.

First, the frequency band of the signal needs to be determined. The signal from sensor #0 is processed by the fast Fourier transform and meshed at 10 kHz intervals (for example, 100–110 kHz). Since the number of points in each frequency band is consistent, all the values in each frequency band are summed and normalized as the energy value of the frequency band. The frequency bands that reach more than 50% of the peak value are reserved, and the bands are defined by the internal frequency as *f1*, *f2*, ... , *fk* (for example, 100–110 kHz is 105 kHz). The total number of the reserved frequency bands is *k*. The weight of the peak frequency band is 1, and the corresponding weights determine the frequency weighting matrix *mf1*, *mf2*, ... , *mfk* of the different frequency bands.

The signals of sensors #0–7 are narrow-band filtered on the above frequency bands to obtain a processed signal matrix. The processed signal is represented as *fk*(*t*,*n*). For leakage excitation in the thin plate, a large out-of-plane displacement is generated [19]. For the out-of-plane motion, compared to the A0 mode, the energy of the S0 mode in the Lamb wave is relatively small [17]. Because the beamforming method mainly focuses on energy, not the arrival time in the TDOA, the wave velocity of the A0 mode is selected as the wave velocity for each bands. The wave velocity matrixes *sf1*, *sf2*, ... , *sfk* at the above frequencies are determined by the dispersion curve. Therefore, the signals of each sensor are delayed and superimposed to obtain the following delayed sum signal *h*(*t*,θ*'*) (corresponds to *g*(*t*,θ*'*) in Section 2.1.):

$$h(t, \theta') = \sum\_{n=1}^{7} \sum\_{k=1}^{K} m f k \times \left[ f\_k(t - \frac{d\_n}{sfk}, n) + f\_k(t, 0) \right],\tag{7}$$

**Figure 2.** Algorithm flowcharts: (**a**) traditional beamforming algorithm and (**b**) frequency weighting matrix beamforming (FWMB) algorithm.

Through the above changes, the FWMB algorithm can improve the traditional beamforming algorithm's ability to locate thin plate leaks and its adaptability to different leaks.

## **3. Experimental Validation**

#### *3.1. Experimental Setup*

The performance of the proposed method was evaluated experimentally. The experimental setup is first introduced in this section. The schematic diagram and photo of the experimental setup are shown in Figures 3 and 4, respectively.

A 5A06 aluminum alloy plate, which is consistent with the material of spacecraft bulkheads, was chosen. The test plate had a diameter of 100 × 100 cm and a thickness of 2 mm. A 1 mm diameter hole was drilled in the center of the board, and this was used to simulate the damage on the spacecraft. Grids with dimensions of 5 × 5 cm were drawn onto the plate to indicate the position of the sensor array.

**Figure 3.** A schematic diagram of the experimental setup.

**Figure 4.** A photo of the experimental setup.

The vacuum pump was used to provide a stable leakage pressure differential to simulate actual leakage conditions. After putting the vacuum-to-plate adapter under the leak hole of the plate, the pumped was turned on until the adapter was attached tightly to the plate. After adjusting the vacuum pump control valve, the subsequent experimental research was able to be started when the measured value of the vacuum degree was less than 10,000 Pa.

The signal acquisition part of the experiment setup was made up by an acoustic emission sensor array, pre-amplifiers, and an acoustic emission instrument. The sensor used in the experiment was a Nano30 sensor (Physical Acoustics Corp., New Jersey, USA) that was 8 mm in diameter and 7.5 mm in height. The Nano30 sensor has a stable frequency response from 100 to 750 KHz. Eight sensors, numbered from 0 to 7, were fixed by special fixtures as 'L,' as shown in Figure 5. The piezoelectric acoustic emission sensors were mounted perpendicular to the plate. The elastic Lamb waves that were caused by the leakage were acquired by the sensor array, and the out-of-plane displacements were mainly measured [20]. The acoustic emission instrument was a DS2-16A instrument (Soft Island Corp., Beijing, China), which can realize multi-channel synchronous acquisition with a 3 MHz sampling rate. The sensors were connected to the pre-amplifiers (Soft Island Corp., Beijing, China). The signals were amplified by 40 dB through the preamplifiers and collected by the acoustic emission instrument.

**Figure 5.** A photo of the sensor array.

## *3.2. Experimental Process*

A hole located at (0, 0) was drilled into the plate, and it is marked as red circular in Figure 6. Once the set-up was equipped with the air exhaust system, the air evacuated from the hole to the vacuum-to-plate adapter could simulate the leakage. Leak-generated ultrasounds spreading through the plate as guided Lamb waves were acquired by the sensors. The sensor array was placed at position A, of which the actual azimuthal direction of the leak was 9◦, and the signal sampling was started as depicted in Section 3.1. The collection of simultaneous leak signals of 8 sensors was repeated 5 times; each data series lasted 0.5 s and was split into 4 segments. Thus, a total of up to 20 segments of data at each position were gained at position A. A series of experiments was done with our array at 9 positions. The positions are marked with green marks in Figure 6. The corresponding azimuthal direction of the leak at positions B, C, D, E, F, G, H and I was 13◦, 20◦, 24◦, 31◦, 45◦, 57◦, 73◦, and 79◦, respectively. The distance from the sensor array to the leak source was maintained at R = 30 cm during the experiment. Thus, 180 segments of data were acquired in total for one leakage hole. In order to test the leakage of different holes to verify the adaptability of the method in this paper, two plates with 1.5 and 2.0 mm holes were used.

**Figure 6.** The position of sensor array in the experiment.

From the FWMB theory, the location of signal source can be determined by applying two sensor arrays. Thus, in the rest part of the paper, the orientation accuracy of the leak source is discussed instead of the positioning accuracy.

#### **4. Results and Discussions**

In this section, the above experimental data are processed and described. First, by taking the experimental data as an example, the calculation method of the related matrix is demonstrated to introduce the method proposed in this paper. Subsequently, the results of the experiment are counted to evaluate the performance of the FWMB algorithm.

## *4.1. The Frequency Weighting Factors with Di*ff*erent Holes*

One advantage of the FWMB algorithm is that it can be adaptively optimized for different leakage holes. In practical applications, the causes of leaks are abundant, so the shapes and sizes of leaks are various. For different leaks, the spectral characteristics may vary greatly. In order to verify this problem, 1.5 and 2.0 mm standard circular leaks were selected for testing.

Because signals below 100 kHz have strong background noises and signals above 500 kHz are very weak, signals between 100 and 500 kHz were chosen for signal processing. The time domain plot, frequency spectrum, and frequency weighting factors for the signals that were sampled at position A of sensor #0 are illustrated in Figure 7.

**Figure 7.** Time domain plot, frequency spectrum, and frequency weighting factors for different holes. (**a**) Time domain plot of a 1.5 mm hole in 100 ms, (**b**) time domain plot of a 2.0 mm hole in 100 ms, (**c**) frequency spectrum of a 1.5 mm hole, (**d**) frequency spectrum of a 2.0 mm hole, (**e**) frequency weighting factor of a 1.5 mm hole, and (**f**) frequency weighting factor of a 2.0 mm hole.

The amplitudes of the acoustic emission signals that were generated by the two leaks in the time domain were close, and from the waveform point of view, they were both difficult to obtain effective information from and were confusing. After performing FFT on the signal, the spectrum energy of different frequency bands was integrated at 10 kHz intervals and normalized to obtain its weighting coefficient matrix. In the figure, the center frequency of the band is used as the label (for example, 100–110 kHz is represented as 105 kHz). It can be seen from the results that the acoustic emission signal that was generated by the 1.5 mm leak was in a higher frequency band. This is consistent with the physical understanding that smaller gaps result in sharper and higher frequency sounds.

In order to further optimize the algorithm, a band energy threshold was set. The bands with energy smaller than 50% of the peak value were discarded in subsequent processing. It can be seen from the figure that the 1.5 mm leak reserved nine frequency bands of 100–190 kHz, while the 2.0 mm leak reserved eight frequency bands of 100–180 kHz. For the above two leaks, although the frequency bands were similar, the weighting coefficients were significantly different. For holes with more specific shapes, larger differences will occasionally be encountered. The above results illustrate the adaptive ability of the method in this paper to different leaks.

## *4.2. The Speed Matrix*

In the previous section, the frequency weighting matrix was obtained. In the FWMB algorithm, another key parameter is the wave velocity matrix. Its role is to accurately match the wave speed in different frequency bands to solve the wave speed error that is caused by the dispersion effect in the flat plate. From the perspective of beamforming theory, wave speed is a critical parameter. An inaccurate wave speed results in the wrong delay, which leads to the wrong location result.

Due to the dispersion characteristics of Lamb waves, signals with different frequency components have different phase velocities. It was assumed herein that the phase velocities were approximately the same in a narrow band of 10 kHz. The numerical solution of the Rayleigh–Lamb equation of the test plate is shown in Figure 8.

**Figure 8.** Dispersion curves of the A0 (symmetric) and S0 (anti-symmetric) mode waves for test plate.

The data were processed into eight or nine frequency bands by narrow band-pass filtering: 100–110, 110–120, ... , 180–190 kHz. The corresponding phase velocities of different frequency components in the dispersion curve were recorded as a velocity matrix. For the low frequency thickness product, between 0.3 and 0.6 MHz · mm, there were only two wave modes (the A0 and S0 modes). Since the A0 mode is dominant in the energy of the wave [17], the A0 mode signal is discussed in the rest of this article. Table 1 lists the corresponding phase speeds *cp* and the weighting factors in different frequency bands.


**Table 1.** Phase speeds of different frequency bands of the A0 mode signals.

Once the phase speeds and the geometry of a sensor array are given, the time delay matrix can be obtained via Equations (1) and (2). With the process in Figure 2, the FWMB algorithm can be produced. When the power gets its maximum, the corresponding assumed angle is the best estimate of the actual direction of arrival.

## *4.3. Location Results*

This section gives a preliminary introduction and discussion of the location results. The 20◦ leakage at sensor array position C for 1.5 mm is taken as an example. Figure 9 shows the results of the traditional beamforming algorithm and the FWMB algorithm proposed in this paper.

**Figure 9.** The results for the 20◦ with different algorithms. (**a**) Traditional beamforming with a 120–130 kHz filter, (**b**) traditional beamforming with a 150-160 kHz filter, and (**c**) result with the FWMB algorithm.

For the traditional beamforming, a separate frequency band is usually used to approximate the wave speed information. In order to have more contrast with the FWMB algorithm, a narrow band similar to that in the FWMB algorithm was selected. Figure 9a shows the results in the case of the 120–130 kHz filtering. It can be seen that a maximum value could be obtained near 20◦, and a good result was obtained. However, the obtained result was 19◦, which, although very close to 20◦, still had some errors. In Figure 9b, another frequency band is used as an example. In this example, two peaks appeared, although there was also a significant peak at 20◦. However, a higher peak appeared at around 70◦ and misjudgment occurred. The reason for the above phenomenon is probably that the traditional beamforming only used one single narrow frequency band, and the energy carried by this frequency band was limited and not necessarily the largest. Therefore, there is a certain possibility that the superposition energy obtained from multiple angles was approximated, thereby obtaining multiple peaks. At the same time, the reflection was also a main interference factor for location. Though the energy was reduced for reflection, it was still a source of interference for positioning. The FWMB method can remove the effect of false peaks by superimposing more frequency band results. The above results show that the traditional beamforming algorithm is prone to errors because it uses less frequency band information and is less robust.

To overcome the above problems, under the conditions of this paper, 64 sets of data overlay were used in the FWMB algorithm to improve the eight sets of data overlay in the traditional beamforming. Therefore, it can be seen in Figure 9c that the result graph of the FWMB algorithm has more fluctuations compared to traditional beamforming. Though the data are not smooth, more accurate results were obtained due to the weighted superposition of the multiple sets of frequency band data. In view of the above preliminary conclusions, we show the data of several other points as further displays.

The results for the leakages in positions I and B are shown in Figure 10. The graphics use a circular coordinate system for easier display and understanding. The red curve and angle value represent the FWMB result, and the blue curve and angle value represent the traditional beamforming result. The frequency band in the label indicates that it is used in traditional beamforming. For the FWMB algorithm, the frequency band is adaptive and does not require manual selection.

**Figure 10.** The results for the leakages in position I (real angle = 79◦) with a 1.5 mm hole and position B (real angle = 3◦) with a 2.0 mm hole by different algorithms in polar coordinates. The polar axis value represents the normalized energy. Red and blue curves represent the normalized energy results with the FWMB algorithm and the traditional beamforming, respectively. (**a**) "I," 1.5 mm, 100–110 kHz; (**b**) "I," 1.5 mm, 110–120 kHz; (**c**) "I," 1.5 mm, 120–130 kHz; (**d**) "B," 2.0 mm, 100–110 kHz; (**e**) "B," 2.0 mm, 110–120 kHz; and (**f**) "B," 2.0 mm, 120–130 kHz.

It can be further seen from the results that the FWMB algorithm optimized the results from two perspectives for traditional beamforming. First, due to the matrix-type multi-data superposition, the discrimination of the results was further enhanced and small errors were resolved, such as within 2◦ in Figure 10a,d,e,f. In addition, through the selection of the energy band and the calculation of the weighting coefficient, the situations where multiple peaks appeared and misjudgment occurred, such as in Figure 10b,c, were improved.

It can also be seen that the FWMB algorithm results also had obvious fluctuations. For example, there is a multiple peak at about 13◦ in Figure 10a–c. This shows that the FWMB algorithm may have misjudged the location, especially when the leakage was near the edge, which was mainly related to the reflection of the boundary. However, due to the large number of FWMB overlays, this problem was compensated for by averaging, and correct results were almost obtained.

## *4.4. Results Statistics and Evaluation*

This section summarizes the results and evaluates the performance of the FWMB algorithm. The statistics of the leakage angle location results for two different holes are shown in Figure 11. The upper and lower ends of the error bar in the figure indicate the maximum and minimum values in the measurement data.

**Figure 11.** Result statistics for different sensor positions (A–I). (**a**) The results statistics of a 1.5 mm hole, and (**b**) the results statistics of a 2.0 mm hole.

As can be seen from the figure, the FWMB algorithm was able to get a close location result for different angles. Judging from the average error, reading errors of 1.33◦ and 1.93◦ were obtained in the holes of 1.5 and 2.0 mm, respectively, and the comprehensive error was 1.63◦. The above results prove that the FWMB method proposed in this paper is adaptive to different leaks. The above advantages are obtained by weighting the energy ratios of different frequency bands in the FWMB method. Due to the use of narrow frequency band, the dispersion problem of Lamb waves is solved, and the wave velocity parameter in the algorithm is more accurate.

#### *4.5. Discussion and Prospection*

The above experiments confirm the improvement of traditional beamforming by the FWMB algorithm. However, the FWMB algorithm still has some issues that need to be addressed in future research.

• The leak hole that was used in this paper was a circular hole. In subsequent research, irregular holes that are caused by real collisions will be tried to more comprehensively evaluate the effectiveness of the method.


In summary, for continuous leakage, the FWMB algorithm provides significant improvements. In subsequent developments, optimization and further testing are needed based on actual applications.

## **5. Conclusions**

This paper proposes a frequency weighting matrix beamforming algorithm to solve the problem of locating continuous leaking sound sources in plate structures. Compared with the traditional beamforming algorithm, this algorithm uses a narrow frequency band to match the wave speed to solve the influence of the Lamb wave dispersion of the thin plate on the signal positioning. Through the energy-weighted method, the algorithm can adapt to different leaks, which has obvious advantages in application. The above scheme was verified through experiments, and 1.5 and 2 mm leaks were tested. The results showed that the frequency weighting matrix beamforming algorithm that is proposed in this paper can obtain a location accuracy of 1.63◦. The above results can be used for the leak detection of on-orbit spacecraft, because the shape of the leak holes that are generated by space debris is inconsistent, so the adaptability of this method can solve the above problems.

**Author Contributions:** L.Q., X.R. and L.S. established the theoretical method; X.L. and L.W. designed the experiments; G.Y. and T.L. performed the experiment; L.Q. and Y.Z. analyzed the data. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Nature Science Fund of China (No. 61903040, 61973227), State key laboratory of precision measuring technology and instruments (Pilab1706), Tianjin Key R&D Program (No. 19YFSLQY00080).

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

## **References**


© 2020 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* **Substructuring of a Petrol Engine: Dynamic Characterization and Experimental Validation**

**Enrico Armentani 1, Venanzio Giannella 2,\*, Roberto Citarella 2, Antonio Parente <sup>3</sup> and Mauro Pirelli <sup>3</sup>**


Received: 22 October 2019; Accepted: 14 November 2019; Published: 19 November 2019

**Abstract:** In this work, the vibration behavior of a 4-cylinder, 4-stroke, petrol engine was simulated by leveraging on the Finite Element Method (FEM). A reduced modelling strategy based on the component mode synthesis (CMS) was adopted to reduce the size of the full FEM model of the engine. Frequency response function (FRF) analyses were used to identify the resonant frequencies and corresponding modes of the different FEM models, and the obtained results were compared with experimental data to get the model validation. Subsequently, modal-based frequency forced response analyses were performed to consider the loads acting during the real operating conditions of the engine. Finally, the impact on vibrations at the mounts, produced by an additional bracket connecting the engine block and gearbox, was also investigated. Both the full and reduced FEM model demonstrated and reproduced with high accuracy the vibration response at the engine mounts, providing a satisfactory agreement with the vibrations measured experimentally. The reduced modelling strategy required significantly shorter runtimes, which decreased from 24 h for the full FEM model to nearly 2 h for the reduced model.

**Keywords:** FEM; component mode synthesis; petrol engine; NVH; FRF

#### **1. Introduction**

Due to the ever-increasing complexity of problem solving and the limitations of power computing, several methods have been developed, since the early days of numerical computation, to improve the efficiency of numerical simulations. Many techniques have been investigated to reduce the computational effort required to solve large static and dynamic problems. Modal analysis techniques have been developed to decouple the large set of ordinary differential equations and to minimize the effort required to solve an iteration of the equations of motion [1,2]. Substructuring methods have been developed to approximate complex big structures as a collection of smaller single components, allowing the simulation process to be performed in a sequence of intermediate computationally lighter steps. These advances proved beneficial, especially in the simulation of dynamic systems, which require the solution of equations of motion over large time intervals and with many individual time steps.

Accurately assessing the dynamic behavior of these elements requires the use of large finite element (FE) models [3–7] to represent the geometry in considerable detail. Assembling the individual sub-components to build up a global FE model of the entire structure results in very large models having a degrees of freedom (DoFs) number which easily exceeds the limits of computer capacity, at least for reasonable runtimes. The question arises whether such FE models can be reduced in

size, preserving at the same time the capability to represent the dynamic characteristics of the entire structure with sufficient accuracy. The substructuring method leverages on dividing a large model into subcomponents that are separately analyzed and afterward re-assembled in a global model, through coupling of their mathematical description.

In recent decades, a variety of methods aimed at a model order reduction of dynamic problems have been developed within the area of structural mechanics, with mode-based methods being the most frequently used. Fairly recently, methods originating from control theory have been employed within structural mechanics. In contrast to mode-based methods, which have an explicit physical interpretation, such modern reduction methods are developed from a purely mathematical point of view, see, e.g., [8–11].

Due to several contributors, such as unbalanced reciprocating and rotating parts, cyclic variations in gas pressure generated by the combustion process, misfire, and inertia forces of the reciprocating parts, significant levels of vibration are induced in internal combustion engines (ICEs). The vibration signals are categorized as torsional, longitudinal, and mixed vibrations. Torsional vibrations are mainly caused by the exertion of cyclic combustion forces within the cylinder as well as the inertial forces of rotating parts, such as the crankshaft, camshaft, and connecting rods. The primary sources of longitudinal vibrations are the unbalanced forces acting on reciprocating and rotating components of the engine which propagate in the structure. The interactions of both longitudinal and torsional vibrations are called mixed vibrations. The inertial forces of rotating components, as well as harmonic combustion forces in the cylinders, eccentric rotation of journal bearing support and flywheel, and transient contact dynamics of the cam/follower, are further contributors to the problem of vibrations in ICEs. Noises generated by these vibrations are transmitted through the engine mounts and chassis to backrests, which can affect passenger comfort and the safety of the vehicle. The level of safety and vibration isolation depends on the amplitude, wave shape, and duration of exposure to these noises. In the review paper [12], different types of vibration in ICEs and their fundamental sources can be found.

In the proper use of the engine mounting system (EMS), when a high transmission ratio is selected, the vibration and noise can be remarkably isolated before they are transferred to the chassis and body of the vehicle. EMS plays an inevitable role in isolating the driver and passenger from the noise, vibration and harshness (NVH) produced by power-intensive engines in modern light-weight vehicles. Therefore, the performance and reliability of these systems require further improvements for a better NVH refinement.

To minimize the transmitting vibrations, the stiffness of the engine block can be enhanced by applying structural modifications such as thickening the crankcase walls, front gear cover, back flywheel cover, as well as adding ribs to its driveshaft.

In a numerical work by Junhong and Jun [13], finite element analysis (FEA) and multi-body analysis (MBA) tools were implemented for modeling the dynamics of the acoustical components of a six-cylinder in-line diesel engine, and also to understand the interaction between excitation mechanisms and noise transmission in the system. To minimize the vibration transferred through the internal paths and to enhance the structural stiffness of the engine, a modification to the design was proposed, namely adding ribs to the drive shaft of the engine.

In the present study, a substructuring technique is applied to a car petrol engine, modelled considering explicitly all its sub-components, i.e., gearbox, exhaust system, alternator, etc.

The objective of the analyses carried out in the present investigation was to validate the adopted approach, based on the component mode synthesis (CMS), to assess the vibrational behavior of an in-line 4-cylinder, 4-stroke, petrol engine, with manual transmission and a total displacement of 1200 cc. A reference solution for the engine was provided by a full FEM model and experimental results from dynamic bench tests were also available as a benchmark.

A frequency response function (FRF) analysis was used to identify the resonant frequencies and mode shapes of the different FEM models. The results were compared in terms of the vibrational response at the mounts, i.e., the gear mount and the engine mount. Experimental results were also available and used for the FEM model validation.

The commercial FEM code MSC Nastran [14] was selected as the FEM solver whereas the commercial code Siemens LMS Virtual Lab [15] was used for the dynamic analyses.

The proposed CMS-based strategy provides the opportunity to circumvent non-disclosure limitations when dealing with original equipment manufacturers (OEMs), since only the subset of few strictly needed data can be exchanged with OEMs without affecting the accuracy of their calculations for those parts of the powertrain designed and manufactured in outsource. Moreover, a further validation concerning the accuracy of such an approach when dealing with complex problems was provided in this work.

## **2. Full FEM Engine Modelling**

The FEM model of the entire engine created by means of the commercial code Altair Hypermesh [16] is shown in Figure 1. Tetrahedral quadratic elements were used to model components such as the engine crankcase, cylinder head, crankshaft, gearbox, etc., whereas quadrilateral quadratic surface elements were used to model thin elements such as the exhaust pipe or the oil pan. Bar elements were used to model the bolts that connect the various components. The average element size was set to 3–6 mm. The final full FEM model comprised nearly 5.7 <sup>×</sup> 106 elements and 4.7 <sup>×</sup> <sup>10</sup><sup>6</sup> nodes.

**Figure 1.** FEM model of the engine with all its sub-components.

In order to validate the full FEM model, two frequency response function (FRF) analyses were performed: one for the vertical bending and one for the lateral bending. An FRF analysis consists of the application of a unitary force in a given node, followed by measurements of the related accelerations in different nodes of the model. Such relevant points selected for the FRF analyses of the engine are shown in Figure 2 for both the vertical and the lateral bending load case. These numerical analyses replicated the experimental tests previously performed and the results are compared in the following section.

**Figure 2.** Full FEM model with details of the relevant points for the FRF analyses for load cases of (**a**) vertical bending and (**b**) lateral bending.

The FEM model was imported in the MSC Nastran to calculate the modal basis of the engine up to a frequency of 1500 Hz (the Lanczos algorithm was used for such a purpose). A damping coefficient equal to 0.03 was considered for all the structural elements. Subsequently, the modal basis was imported in the Siemens LMS Virtual Lab environment, with input and output points defined for either the vertical (Figure 2a) or lateral (Figure 2b) bending. Finally, FRF and modal-based frequency response analyses were performed to calculate the accelerations at the relevant points.

A free-free condition was simulated in all the presented cases. As a matter of fact, the engine mounts are capable of decoupling the dynamic engine behavior from the supporting structure in the whole frequency range considered.

Moreover, the impact of the crankcase–gearbox bracket shown in Figure 3 on the acceleration levels at the previously mentioned relevant points was assessed. Such evaluation was needed to understand whether the impact on the dynamic behavior of the engine was necessary, since its removal would have a positive impact on cost reduction.

**Figure 3.** Close up of the crankcase–gearbox bracket under analysis.

## **3. FRF Results of the Full FEM Model**

The accelerations were calculated at the engine and gear mounts for the two analyses of vertical and lateral bending (Figure 2). The ratio between acceleration and input force is shown in Figures 4 and 5. In particular, Figure 4 shows results for the vertical bending at both mounts, with and without the engine bracket shown in Figure 3. It is worth noting that the impact of the engine bracket was much more relevant for the vertical bending load case (Figure 4), whereas it played a minor role for the lateral bending load case (Figure 5). This result was expected due to the geometry of the bracket, which allows it to stiffen the engine–gearbox connection, especially when undergoing vertical loads.

**Figure 4.** *Cont*.

**Figure 4.** Acceleration/force ratio for the full FEM model for the vertical bending load case with and without the crankcase–gearbox bracket (shown in Figure 3) at: (**a**) engine mount; (**b**) gear mount.

**Figure 5.** Acceleration/force ratio for the full FEM model for the lateral bending load case with and without the crankcase–gearbox bracket (shown in Figure 3) at: (**a**) engine mount; (**b**) gear mount.

From Figure 4, it is possible to see that the bracket introduction is effective in increasing the first natural frequency, in such a way as to avoid any activation of such mode in the whole engine operating range, up to 6000 rpm. As a matter of fact, the second order engine's highest frequency is equal to 200 Hz, well below the first natural mode whose frequency is now nearly 220 Hz.

From Figure 5, it is possible to see that even if the increase in the first natural frequency is not sufficient to overcome 200 Hz, the corresponding peak is severely smoothed in such a way that it no longer represents a problem, whereas the second more relevant peak is relegated to a frequency much higher than 200 Hz.

A comparison between the numerical and experimental results is also provided in Figure 6. The FRF results show a satisfactory correlation between experimental outcomes in terms of peak frequencies. The discrepancy in terms of magnitudes was judged as physiological for this kind of analyses because of uncertainties on the correct damping value to be used for the simulation, and because, when modal reduction techniques are adopted it is difficult to make the experimental and numerical exciting/accelerometer points coincident. Figure 6a shows a discrepancy between the peak frequencies of nearly 5 Hz for the model with the engine bracket, whereas Figure 6b shows a discrepancy of nearly 10 Hz between the peak frequencies for the model without the engine bracket. In both cases, such approximations are judged acceptable.

**Figure 6.** Numerical/experimental comparison on the acceleration/force ratio measured at the engine mount for the vertical bending load case: (**a**) with and (**b**) without the crankcase–gearbox bracket shown in Figure 3.

## **4. Frequency Response Results of the Full FEM Model**

Modal-based frequency response analyses were performed for the full FEM engine model, considering the real loads occurring during the operation of the engine. In particular, the wide open throttle (WOT) condition was simulated, thus the maximum intake of air–fuel mixture, occurring when the throttle is completely opened, was considered as the engine working condition. Moreover, vertical loads representing the pistons and roll loads (from crankshaft rotation) were considered. For the former, only loads corresponding to the second engine order were needed (being the remaining negligible), whereas the fourth and sixth order were also added for the latter. Such loads were applied on a RBE3 element that in turn distributed the load onto the elements of interest. The objective of these analyses was to simulate the impact of the aforementioned bracket on the vibrations measured at the mounts during the engine operation. Experimental measurements of the mount vibrations were also available and were used here for validation.

Figure 7 shows the magnitude of the numerical displacement evaluated at the engine and gear mounts, with and without consideration of the presence of the bracket. It is worth noting that the adoption of the bracket has a positive effect at the engine mount just inside the range 5000–5500 rpm, whereas it seems to be ineffective elsewhere. On the contrary, a sensible positive bracket impact on the gear mount total displacement is evident above 4500 rpm. Apparently, only the bracket turns out to increase the magnitude of gear mount vibrations in the 5500–6000 rpm range because, as shown in the following, the model prediction in such a range, with reference to the gear mount, is inaccurate. Considering the experimental outcomes, as expected, the impact of the bracket introduction is also positive in the range 5500–6000 rpm: it is therefore sufficient to compare the blue testing lines of Figure 8a,b to understand that the total displacement is always lower when considering the bracket addition.

**Figure 7.** Displacement magnitude calculated at the (**a**) engine and (**b**) gear mounts considering the realistic operating loading condition.

Figures 8 and 9 show the experimental/numerical comparisons in terms of the displacements at the gear and engine mounts respectively. The setup for the experimental modal analyses is shown in Figure 10, which also highlights the supporting structure adopted to hold the engine.

As previously anticipated, an acceptable correlation was obtained between the data but not for the gear mounts without the bracket (Figure 8b). A possible explanation could be related to the approximations inherent in the modelling of interface connections between the engine parts, where nonlinear contact conditions are skipped in order to reduce the computational burden but with the consequence of possible unrealistic interpenetration between interface surfaces. This aspect can become critical at high frequency (or namely at high engine rpm) as shown in Figure 11, where it is possible to observe that the gear mount presented excessive relative motions with reference to the gearbox at 210 Hz (referring to the second order, this corresponds to engine rpm slightly higher than 6000). Some remaining discrepancies can be ascribed to the uncertain damping attributed to the structural elements and to the simplified modelling of the bolts that interconnect the different parts of the engine. In any case, the full FEM modelling was considered as validated.

**Figure 8.** Experimental/numerical comparison in terms of displacements at the gearbox mount (**a**) with and (**b**) without consideration of the bracket highlighted in Figure 3.

**Figure 9.** Experimental/numerical comparison in terms of displacements at the engine mount (**a**) with and (**b**) without consideration of the bracket (highlighted in Figure 3).

**Figure 10.** Experimental setup for modal analyses with highlights of the supporting structure.

**Figure 11.** Evidence of the limitations given by the modelling strategy of the contacts, with reference to the relative motion between gear mount and gearbox at 210 Hz.

#### **5. Model Reduction**

The reduced FEM model was created starting from the full FEM engine model previously presented and validated. This was achieved by means of the ANSA MetaPost [17] code together with the pre-processor Altair Hypermesh and MSC Nastran.

The reduction of the model was performed for the whole engine apart from the sub-components. Namely, the reduced model comprised a fully reduced engine model with the addition of the engine's non-reduced sub-components (still fully modelled by FEM). This modelling strategy was preferred since it allowed efficient accommodation of the OEMs (Original Equipment Manufacturers) requirements, as they generally require the explicit modelling of their specific sub-components when assembled to an engine. From this standpoint, the model reduction of the engine allowed the significant reduction of the size of the FEM model, and also enabled the possibility of sharing the model without disclosing proprietary information about the engine design.

Sub-components were removed from the full FEM model and PLOTEL elements (PLOT ELements) [14] were used to display the engine shape and size. PLOTEL elements are MSC Nastran dummy elements used only for display purposes and without significance from the dynamic standpoint. Such elements were adopted to represent the correct shape in the relevant positions, e.g., load application points or bolt positions, and also to visualize the correct engine shape and size.

The total mass and stiffness of the engine was associated with an RBE3 [14] element positioned in the cylinder head, whereas further RBE3 elements were also introduced to link the sub-components to the engine, in correspondence with the connecting bolts (these elements allowed transmission of displacements at the relevant positions where sub-components are directly connected).

The final reduced model is shown in Figure 12.

**Figure 12.** (**a**) Engine model built up with PLOTEL elements; (**b**) highlight of RBE3 elements.

The model built up with PLOTEL elements is shown in Figure 12a. It comprised 1885 elements and 1755 nodes (whereas the full FEM model required nearly 5.7 <sup>×</sup> 106 elements and 4.7 <sup>×</sup> 106 nodes). Such a model was then imported in the ANSA MetaPost code, in which the modal model was built up by considering the modes of the full FEM model as input data for the reduced model. Consequently,

further RBE3 elements were used to link the sub-components at the bolt position, see Figure 12b. The final reduced model comprising all the FEM sub-components is shown in Figure 13.

**Figure 13.** Reduced model of the engine with all sub-components modelled by FEM.

FRF analyses were used to calculate the accelerations at the mounts for the full and reduced FEM model. In this case, unitary forces were applied on the sub-components and the corresponding accelerations were obtained (in the same position of excitation). For instance, Figure 14 shows the points considered in the FRF analyses on the gear mount (Figure 14a), on the alternator (Figure 14b) and on the intake manifold.

**Figure 14.** Relevant points (yellow dots) on the (**a**) gear mount; (**b**) alternator and (**c**) intake manifold considered in the FRF analyses.

## **6. Results of the Reduced FEM Model**

Figure 15 shows the ratios of acceleration over force calculated for the two FRF analyses performed for the alternator and for the gear mount: a quite satisfactory matching was obtained.

**Figure 15.** Comparisons of the acceleration/force ratio in the three directions for the full and reduced models at: (**a**–**c**) gear mount; (**d**–**f**) alternator; (**g**–**i**) intake manifold.

Figure 16 shows three examples of the modal shapes for the full FEM model and the reduced model for three components.

**Figure 16.** Comparisons of the modal shapes for (**a**,**c**,**e**) the full FEM and (**b**,**d**,**f**) the reduced FEM model for: (**a**,**b**) gear mount, (**c**,**d**) alternator, (**e**,**f**) intake manifold.

Figure 16 shows an example of the modal shapes for the full FEM model and the reduced model: again, a quite satisfactory matching was obtained.

The reduced FEM modelling allowed the calculation of the accelerations at the sub-components with a reduced computational burden. In particular, the runtime for the full FEM model was equal to nearly 24 h, decreasing to nearly 2 h for the reduced FEM model. The reduced model was then able to accurately replicate the engine's vibrational behavior and consequently could be adopted for this kind of analysis. Moreover, since the sub-components are generally designed by external partners who require the analysis of the sub-component assembled to the engine, the adoption of the model reduction can overcome the key concern of not sharing proprietary design solutions.

## **7. Conclusions**

A 4-cylinder petrol engine was simulated from the dynamic standpoint by leveraging on the finite element method (FEM). In particular, a reduced modelling strategy based on the component mode synthesis (CMS) was adopted to reduce the size of the full FEM model of the engine.

The FEM model of the engine, comprising all its sub-components, was preliminarily characterized from the vibration standpoint; subsequently, the CMS was adopted in order to reduce the FEM model DoFs.

FRF analyses were used to reproduce the vibration response of the engine and the corresponding acceleration levels between the models and bench test data were compared, showing a sound agreement. Moreover, frequency response analyses were also performed to replicate the vibrations of the engine during the operation and a satisfactory comparison was obtained against corresponding bench test data.

The adopted reduced modelling strategy turned out to be effective in lowering the computational burden from 24 h to nearly 2 h, keeping at the same time an accurate replication of the engine vibration behavior.

The reduced FEM model was demonstrated to reproduce with high accuracy the vibration response at the engine mounts, providing a satisfactory agreement with the vibrations measured experimentally and with the outcomes of a full FEM model.

Ways to further improve the modelling of the bolts and of the interface contacts between the engine parts are still under investigation.

**Author Contributions:** A.P., M.P. and E.A. conceived and performed the simulations; V.G. wrote the paper, R.C. supervised the work.

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

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

## **References**


© 2019 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/).
