# Attenuation of Line Spectra and Vibration Based on a Chaotification Method and Global Analysis

Authors: Zhao-Lun Zuo; Xiang Yu; Shi-Jian Zhu
Journal: Shock and Vibration
 (2022)

---

## Abstract

The chaotification method can be used to weaken the line spectra component of the ship's radiated noise. The main obstacle that plagues the engineering realization of this innovative concept is how to realize the small-amplitude chaotic motion of the nonlinear vibration isolation system (NVIS) in a large parameter range and how to maintain the chaotic state under variable operating conditions. In this paper, the dimensionless dynamic equation of the NVIS with a flexible base is derived. A control method similar to generalized synchronization is used to realize the continuous chaos of the system. The maximum Lyapunov exponent, the conditional Lyapunov exponent, and disturbances are introduced to verify the effectiveness of the control method and its stability and robustness in a large range of parameters. Aiming at the coexistence of multiple attractors, a new global analysis method is proposed, which is used to study the distribution of basin of attraction (BA). According to the changes of BA before and after the system is controlled, the BA is partitioned. Initial conditions from different partitions are selected to study the attenuation effect of line spectra and vibration, and the favorable initial conditions for realizing the small-amplitude chaotic motion of the system are given. The method may be extended to the global analysis of other controlled systems.


---

## Body

## 1. Introduction

The chaotification (chaos anticontrol) method [1–4] has many applications in engineering, one of which is to weaken the line spectra component of the ship's radiated noise. This method [5] uses an NVIS in a chaotic state to convert the line spectra input into a broadband continuous spectra output, so that the energy originally concentrated at the line spectra is dispersed to the entire frequency band. The main obstacle that plagues the engineering realization of this innovative concept is how to realize the small-amplitude chaotic motion of the NVIS in a large parameter range and how to maintain the chaotic state under variable operating conditions.The chaotic synchronization method can convert the system from the periodic state to the chaotic state [6, 7]. Many scholars apply this method to the line spectra attenuation of marine mechanical NVIS. Yu et al. [8], based on the generalized chaotic synchronization principle, used an external chaotic system to drive some parameters of the NVIS, thus making the system always maintain chaotic state. Since the response spectra of chaotic system under harmonic excitation is continuous and weakened, the line spectra is effectively reduced. Li et al. [9] used the Lorenz system family to chaotify the NVIS and proposed a strategy to obtain the optimal control gain to ensure the quality of chaotification. Due to the inherent characteristics of chaotic behavior, the amplitude of chaotic vibration is usually greater than that of its specified periodic motion. In order to overcome the conflict between line spectra attenuation based on chaotic vibration response and vibration isolation capability, Wen et al. [10] developed a method based on modified projective synchronization to improve vibration isolation performance and utilized broadband frequency chaotic signal to achieve line spectra reduction.However, as a typical nonlinear system, chaotic synchronization systems often present complex response characteristics, including multivalue mapping relationships, namely, multiattractor coexistence [11]. Due to the jumping phenomenon, the large-amplitude and small-amplitude attractors that may appear in a mechanical system depend on the initial conditions. When the system operates on a large-amplitude attractor, namely, the resonant branch, it is very unfavorable to the engineering implementation of various chaotification methods. Therefore, it is necessary to predict and analyze the global behavior of the system and determine the initial conditions that generate different response behaviors, including (1) determining the type and amplitude of each attractor and (2) determining the size and distribution of the BA associated with each attractor. Through global analysis and drawing the diagram of the BA, the initial conditions for generating various response behaviors can be visually displayed. The most intuitive method to determine the BA is the Poincaré mapping method [12, 13], which discretizes the time domain of the system for sampling and preserves the continuity of the phase space. The advantage of this method is the high accuracy in determining the distribution of the BA, but it has the limitation of large computational effort. The simple cell mapping method proposed by Hsu et al. [14] discretizes the time domain and phase space to ensure a certain accuracy and at the same time improve the computational efficiency. Subsequently, various cell mapping methods such as interpolated cell mapping [15] and generalized cell mapping [16] have been developed and applied to the dynamic analysis of different systems. These cell mapping methods need to determine the mapping sequence of each cell in the N-dimensional cell state space, which leads to an exponential increase in the computational effort as the system dimension increases. Therefore, for high-dimensional systems, the computational efficiency of the various cell mapping methods is not as good as that of the Poincaré mapping method.In this paper, a control method similar to generalized synchronization [8] is used to realize the chaotification of the controlled NVIS, and its effectiveness, stability, and robustness are verified. Based on path tracking, phase trajectory reconstruction, and constrained state space selection, a new global analysis method is proposed to study the distribution of BA of the system. Compared with the point mapping method and the cell mapping method, this method splits the process of determining attractor and BA, reduces the cycle operation, and greatly improves the computational efficiency. The types of coexisting attractors and their distribution in phase space are visualized by path tracking and phase trajectory reconstruction. By choosing the analysis plane and the restricted state space, the calculation is independent of the dimensionality of the system, which is suitable for high-dimensional systems. According to the changes of BA before and after the system is controlled, the BA is partitioned. The initial conditions from different partitions are selected to study the attenuation effect of line spectra and vibration, and the favorable initial conditions for realizing the small-amplitude chaotic motion of the system are given.The paper consists of five sections. Following is an introduction of the background in Section1. The dimensionless dynamic equation of the NVIS on a flexible base is derived in Section 2. In Section 3, a control scheme of the unidirectional coupling mode based on generalized synchronization is used to realize the persistent chaotic motion of the system; the stability and robustness of the scheme are demonstrated. In Section 4, the process of global analysis is introduced; the types of coexisting attractors and their distribution ranges in phase space are visualized by path tracing and reconstruction of the phase trajectory. According to the changes of BA before and after the system is controlled, the BA is partitioned. The initial conditions from different partitions are selected to study the attenuation effect of line spectra and vibration, and the favorable initial conditions for realizing the small-amplitude chaotic motion of the system are given. In Section 5, some conclusions are drawn with further works discussed.


## 2. Dynamical Model

In practical applications, it is necessary to consider the flexibility of the base of the marine mechanical NVIS. When only the first-order mode is considered, the flexible base can be simplified as a lumped mass, and the entire NVIS can be regarded as a two-degree-of-freedom (DOF) mass-spring system. Schematic diagram of the controlled NVIS is shown in Figure1; M1 and M2 represent the device and the base, respectively. M1 is supported by a vibration isolator that combines a linear damper and a nonlinear spring with cubic nonlinearity. M2 uses linear dampers and linear springs that are much harder than nonlinear springs to connect to the fixed plane. Install the actuator between the device and the base. The physical picture of the control system of the double-layer vibration isolation test bench is shown in Figure 2.Figure 1
Schematic diagram of the controlled NVIS with flexible base.Figure 2
Control system of the double-layer vibration isolation test bench.When the spring is in the state of natural length, the position of the mass is selected as the origin of the coordinate, and the dynamic equation of the system can be given as follows:(1)M1Z¨1=−C1Z˙1−Z˙2−K1Z1−Z2−K3Z1−Z23+F0cosΩT−M1g,M2Z¨2=−C2Z˙2−K2Z2+C1Z˙1−Z˙2+K1Z1−Z2+K3Z1−Z23−M2g,where M1 is the mass of the device; M2 is the mass of the base; K1, K3, and C1 are linear stiffness, cubic nonlinear stiffness, and damping coefficient of nonlinear isolator, respectively. K2 and C2 are stiffness coefficients of linear spring and damping coefficients of damper between base and fixed plane, respectively. Z1 and Z2 are displacement variables of device and base, respectively.When the spring is compressed and in static equilibrium, the position of the mass is chosen as the origin of the new coordinates, and the relationship between the new and old coordinates is(2)Z1=Y1−h1Z2=Y2−h2,where Y1 and Y2 are the new displacement variables of the device and the base; h1 and h2 are the compressions of nonlinear spring and linear spring. The following force relationship can be used to eliminate the gravity term at the right end of (1)(3)M1g=K1H+K3H3M1g+M2g=K2h2,where H=h1−h2. Substituting (2) and (3) into (1), we get(4)M1Y¨1+C1Y˙1-Y˙2+K1+3K3H2Y1-Y2-3K3HY1-Y22+K3Y1-Y23=F0cosΩT,M2Y¨2+C2Y˙2+K2Y2=C1Y˙1-Y˙2+K1+3K3H2Y1-Y2-3K3HY1-Y22+K3Y1-Y23,let K0=K1+3K3H2 ; introduce new time scale T1=M1/K0, length scale L1=K0/K3, and dimensionless parameters yi=Yi/L1, t=T/T1, i=1,2; we get(5)dYidT=L1T1⋅dyidt,d2YidT2=L1T12⋅d2yidt2,i=1,2.Substituting (5) into (4), the first-order form of the dimensionless dynamic equation can be obtained(6)y˙1=y2,y˙2=−ξ1y2−y4−k1y1−y3+λy1−y32−y1−y33+fcosωt,y˙3=y4,y˙4=−εξ2y4−εk2y3+εξ1y2−y4+εy1−y3−ελy1−y32+εy1−y33,where ξi=Ci/MK0, ε=M1/M2, λ=3H/L1, k1=1, k2=K2/K0, f=F0/K0L1, ω=ΩT1, i=1,2. Set the parameters of the system as follows: ε=0.2, λ=1, k2=2, ω=3.9311 , ξi=0.8 , i=1,2.


## 3. Control Principle and Validity

### 3.1. Control Principle

In generalized synchronization systems, the relationship between subsystems can be divided into bidirectional coupling [17, 18] and unidirectional coupling [19, 20]. Bidirectional coupling means that the drive and response systems influence each other. The unidirectionally coupled system is also known as the master-slave system; its essence is to adjust the response system according to the output of the drive system. In this mode, the output of the response system complies with the output of the drive system in some form without affecting the drive system.In this paper, we adopt the unidirectional coupling mode and consider the following general form:(7)x˙=Fx,y=Gy,kPx,where x=x1,x2,…,xn and y=y1,y2,…,ym are the state vectors in the n-dimensional space Rn and m-dimensional space Rm, respectively; F:Rn⟶Rn and G:Rm⟶Rm are the vector fields representing the drive system and the response system, respectively. P:Rn⟶Rm denotes the coupling vector field between the subsystems, and k is the scalar coefficients that represent the coupling strength.The NVIS with a flexible base is used as the response system and the drive system is selected as the following single DOF harmonic excitation oscillator:(8)x˙1x˙2=x2−ux2−wx1−vx12−x13+dcosω1t.The parameters are set as follows:ω1=3.9311, w=4, v=3, u=0.15, d=9. Without loss of generality, the drive system’s initial conditions (x1, x2) are set to (0, 0). As shown in Figure 3(a), the maximum Lyapunov exponent of the drive system is greater than 0, indicating that the system is chaotic. Taking x2 as the chaotic driving signal, its time history is shown in Figure 3(b).Figure 3
(a) Lyapunov exponents of the drive system; (b) time history of the drive signalx2.
(a)(b)Stiffness and damping can be controlled in engineering. In this paper, the dimensionless linear stiffnessk1 is taken as a variable parameter and the chaotic output x2 is used to perturb it to realize the persistent chaotic state of the NVIS. Consider the simplest case k1 = kx2 and let the coupling coefficient k = 0.1. The dynamical equation of the response system can be rewritten as(9)y˙1=y2,y˙2=−ξ1y2−y4−kx2y1−y3+λy1−y32−y1−y33+fcosωt,y˙3=y4,y˙4=−εξ2y4−εk2y3+εξ1y2−y4+εy1−y3−ελy1−y32+εy1−y33.


### 3.2. Validity, Stability, and Robustness

As shown in Figures4(a) and 4(b), the maximum Lyapunov exponent of the response system is less than 0 before being driven and greater than 0 after being driven, indicating that the system changes from the periodic state to the chaotic state.Figure 4
Lyapunov exponents of the response system (a) before and (b) after being driven.
(a)(b)The stabilization of the synchronization system means that all transverse motions of the synchronization manifold must be suppressed, which can be proved by conditional Lyapunov exponent (CLE). If all CLEs are negative, it will prove that the synchronization manifold is stable; otherwise, it is unstable [21–23]. As shown in Figures 5(a) and 5(b), when k = 0.1, the maximum CLE of the controlled system λmax  = -0.063; when k is in the range of 0–0.5, λmax is less than 0. It can be concluded that the controlled system tends to be stable in a wide range of parameters.Figure 5
CLE andλmax of the controlled systems. (a) CLE where k = 0.1. (b) λmax where k = 0–0.5.
(a)(b)To analyze whether the controlled system is in transient chaos and investigate the robustness of the control, set two groups of different initial conditions (0, −5, 2, 0), (0, −5, 2.5, 0.5), (0, −5, 1.5, −0.5) and (0, −5, 0, 0), (0, −5, 0.5, 0.5), (0, −5, -0.5, 0.5) to simulate the control system subjected to disturbances. In the NVIS of marine machinery, the vibration is transmitted to the hull through the base. Therefore, the following analysis focuses on the displacement responsey3 of the base. As shown in Figures 6(a) and 6(b), the controlled system does not turn into the periodic state after the transient process but remains chaotic, thus proving that the system is not in transient chaos. The controlled system is asymptotically stable to the undisturbed attractor under slight disturbance, indicating that the applied control has certain robustness.Figure 6
Time histories ofy3 under different initial conditions: (a) small amplitude; (b) large amplitude.
(a)(b)


## 3.1. Control Principle

In generalized synchronization systems, the relationship between subsystems can be divided into bidirectional coupling [17, 18] and unidirectional coupling [19, 20]. Bidirectional coupling means that the drive and response systems influence each other. The unidirectionally coupled system is also known as the master-slave system; its essence is to adjust the response system according to the output of the drive system. In this mode, the output of the response system complies with the output of the drive system in some form without affecting the drive system.In this paper, we adopt the unidirectional coupling mode and consider the following general form:(7)x˙=Fx,y=Gy,kPx,where x=x1,x2,…,xn and y=y1,y2,…,ym are the state vectors in the n-dimensional space Rn and m-dimensional space Rm, respectively; F:Rn⟶Rn and G:Rm⟶Rm are the vector fields representing the drive system and the response system, respectively. P:Rn⟶Rm denotes the coupling vector field between the subsystems, and k is the scalar coefficients that represent the coupling strength.The NVIS with a flexible base is used as the response system and the drive system is selected as the following single DOF harmonic excitation oscillator:(8)x˙1x˙2=x2−ux2−wx1−vx12−x13+dcosω1t.The parameters are set as follows:ω1=3.9311, w=4, v=3, u=0.15, d=9. Without loss of generality, the drive system’s initial conditions (x1, x2) are set to (0, 0). As shown in Figure 3(a), the maximum Lyapunov exponent of the drive system is greater than 0, indicating that the system is chaotic. Taking x2 as the chaotic driving signal, its time history is shown in Figure 3(b).Figure 3
(a) Lyapunov exponents of the drive system; (b) time history of the drive signalx2.
(a)(b)Stiffness and damping can be controlled in engineering. In this paper, the dimensionless linear stiffnessk1 is taken as a variable parameter and the chaotic output x2 is used to perturb it to realize the persistent chaotic state of the NVIS. Consider the simplest case k1 = kx2 and let the coupling coefficient k = 0.1. The dynamical equation of the response system can be rewritten as(9)y˙1=y2,y˙2=−ξ1y2−y4−kx2y1−y3+λy1−y32−y1−y33+fcosωt,y˙3=y4,y˙4=−εξ2y4−εk2y3+εξ1y2−y4+εy1−y3−ελy1−y32+εy1−y33.


## 3.2. Validity, Stability, and Robustness

As shown in Figures4(a) and 4(b), the maximum Lyapunov exponent of the response system is less than 0 before being driven and greater than 0 after being driven, indicating that the system changes from the periodic state to the chaotic state.Figure 4
Lyapunov exponents of the response system (a) before and (b) after being driven.
(a)(b)The stabilization of the synchronization system means that all transverse motions of the synchronization manifold must be suppressed, which can be proved by conditional Lyapunov exponent (CLE). If all CLEs are negative, it will prove that the synchronization manifold is stable; otherwise, it is unstable [21–23]. As shown in Figures 5(a) and 5(b), when k = 0.1, the maximum CLE of the controlled system λmax  = -0.063; when k is in the range of 0–0.5, λmax is less than 0. It can be concluded that the controlled system tends to be stable in a wide range of parameters.Figure 5
CLE andλmax of the controlled systems. (a) CLE where k = 0.1. (b) λmax where k = 0–0.5.
(a)(b)To analyze whether the controlled system is in transient chaos and investigate the robustness of the control, set two groups of different initial conditions (0, −5, 2, 0), (0, −5, 2.5, 0.5), (0, −5, 1.5, −0.5) and (0, −5, 0, 0), (0, −5, 0.5, 0.5), (0, −5, -0.5, 0.5) to simulate the control system subjected to disturbances. In the NVIS of marine machinery, the vibration is transmitted to the hull through the base. Therefore, the following analysis focuses on the displacement responsey3 of the base. As shown in Figures 6(a) and 6(b), the controlled system does not turn into the periodic state after the transient process but remains chaotic, thus proving that the system is not in transient chaos. The controlled system is asymptotically stable to the undisturbed attractor under slight disturbance, indicating that the applied control has certain robustness.Figure 6
Time histories ofy3 under different initial conditions: (a) small amplitude; (b) large amplitude.
(a)(b)


## 4. Global Analysis

### 4.1. Process of Global Analysis

In this paper, a new global analysis method is proposed to study the distribution of BA: Tracking the path of the attractor based on the Poincaré mapping method and the pseudo-arc-length continuation method and reconstructing the phase trajectory of different attractors through the point on the path. In this way, the types of coexisting attractors and their distribution in the phase space can be visualized, providing a priori information for the determination of BA. Select the analysis plane and divide uniformly small grids, and determine the final behavior of the center points of all grids based on the information of the coexisting attractors. The method is described in detail below.In general, anN-dimensional continuous-time nonautonomous system can be transformed into an autonomous one by taking time t as a new one-dimensional variable(10)x˙=fx,s,tt˙=1x,s,t∈RN×RL×R+,where xt∈RN is the state variable, s∈RL is the parameter variable, t is the time variable, and f:RN×RL×R+⟶RN is the vector field. If fx,s,t+T=fx,s,t is satisfied, then a global Poincaré section Σ=x,s,t∈RN−1×RL×R+|t=t0 and Poincaré mapping Px0=xx0,t0+T can be defined. Setting arbitrary initial condition x0,t0 that all phase trajectories of the system intersect with the section Σ transversely, the Poincaré mapping has the same time interval τ=T for all points x0∈Σ. If the Poincaré mapping is continuously differentiable, when the variable parameter changes continuously in small steps, the attractor will move smoothly on the section Σ , so that the path of the attractor can be tracked [24].In this paper, the pseudo-arc-length continuation algorithm is used to track the path of the attractor, which is a prediction-correction method with the following basic ideas: First, determine the variable parameters and iterative step length, and give the preselected parameters and the direction of parameter change. Set the initial conditions; different initial conditions correspond to different Poincaré sections. Then start the prediction-correction process, use the current parameters to predict the steady-state trajectory of the system, and use a convergence algorithm such as the Newton-Raphson method to refine the prediction. The path tracking of the attractor is achieved through an iterative process of modifying parameter values. The path of the attractor under different initial conditions is obtained by changing the direction of the parameter change. It is worth noting that path tracking can only qualitatively analyze the attractor types under different initial conditions but cannot determine the amplitude information. This is because the amplitude of the path is closely related to the selected Poincaré section, so the true amplitude of the attractor cannot be displayed. To this end, we use points on the path to reconstruct the phase trajectories of different attractors and then obtain the distribution of attractors in the phase space.Schematic diagram of the global analysis method is shown in Figure7. The distribution of the coexisting attractors in the phase space is used as a basis for judging the different attractors. The analysis plane P:y1,…yi,yi+3…yn=0,yi+1,1≤yi+1≤yi+1,2,yi+2,1≤yi+2≤yi+2,2 and the restricted state space E:yi+1,3≤yi+1≤yi+1,4,yi+2,3≤yi+2≤yi+2,4 are selected to distinguish different attractors. The analysis plane is finely and evenly meshed. The center point x0 of the initial grid is selected as the initial condition, and the phase trajectories of the attractors are obtained by numerical integration methods such as fourth-order Runge-Kutta. According to whether the phase trajectory is in the restricted state space E, the attractor corresponding to the point x0 is determined. After that, a new center point is selected on the analysis plane to carry out the above analysis process until the final behaviors of all center points can be determined; then all the attractors and their BA in phase space are obtained.Figure 7
Schematic diagram of the global analysis method.


### 4.2. Coexisting Attractors

The model of NVIS established in this paper is a typical nonlinear system with cubic nonlinear stiffness. Moreover, the nonlinear system has a jumping phenomenon. In some intervals of the bifurcation parameters, the same bifurcation value corresponds to three different values of the response amplitude. Among them, the intermediate value corresponds to an unstable state, and the large and small amplitudes that may appear in the mechanical system depend on the initial conditions. In order to analyze the internal mechanism of the coexisting attractors in this system, the path of the attractor is tracked by bifurcation analysis, and the continuation direction is changed to determine the parameter range of coexisting attractors.Poincaré sectionΣ=x,s,t∈RN−1×RL×R+|t=0 and Poincaré mapping Px0=xx0,T are chosen such that all intersections of the phase trajectory with Σ have the same time interval τ=T=2π/ω. With the excitation amplitude f as the variable parameter and taking a step size of 0.01, path changes of attractors for forward continuation (the value of f varies from 0 to 50) and backward continuation (the value of f varies from 50 to 0) are shown in Figures 8(a) and 8(b), respectively. It can be seen that the limit cycle of the system flow is converted into a fixed point on the Poincaré section by Poincaré mapping, and there exists a jump interval f∈15.45,24.37 for the fixed point under different initial conditions. In this parameter interval, there exist two coexisting period-1 attractors for the response system.Figure 8
Path changes of attractors for (a) forward continuation and (b) backward continuation.
(a)(b)Without loss of generality, the pointf = 20 is chosen to reconstruct the phase trajectories of the attractors. As shown in Figures 9(a) and 9(b), when the response system is not driven, set the initial conditions as (0,−5,0,−5) and (0,−5,0,0), respectively, the phase trajectories are two closed loops with significantly different amplitudes, and the power spectra density (PSD) has obvious characteristics of line spectra. It can be inferred that there is a large-amplitude period-1 attractor and a small-amplitude period-1 attractor coexisting in the system; the number and type of the attractors verify the results in Figure 8. As shown in Figures 10(a) and 10(b), phase trajectories of coexisting attractors when the response system is driven are entangled and spread across a certain phase space, and PSD is also converted from line spectra to continuous spectra, which verifies the effectiveness of the control.Figure 9
(a) Phase trajectories and (b) PSD of coexisting attractors when the response system is not driven.
(a)(b)Figure 10
(a) Phase trajectories and (b) PSD of coexisting attractors when the response system is driven.
(a)(b)


### 4.3. Partitioning of BA

The analysis planey1=0,y2=0,−10≤y3≤10,−10≤y4≤10 is chosen, and the small-amplitude attractor is confined in the state space E:−0.05≤y3≤0.4,−0.4≤y4≤0.4, which is used as the basis for judging different attractors. On the analysis plane, y3,y4∈−10,10 is divided into 100 equidistant parts, respectively, and 104 grids with the same spacing are obtained. The center point of each grid is used as the initial condition, and the corresponding phase trajectory is obtained by numerical integration. If the phase trajectory is in the state space E, the behavior of the point is determined to be a small-amplitude attractor; otherwise it is a large-amplitude attractor.When the response system is not driven/driven, distributions of BA are shown in Figures11(a) and 11(b), respectively, where the blue and red regions correspond to the large-amplitude attractor and small-amplitude attractor, respectively. It can be seen that BA corresponding to different attractors have changed, such as the points (-1.7, 0.1) and (9.3, 8.5) marked in the figure.Figure 11
Distribution of BA when the response system is (a) not driven or (b) driven.
(a)(b)The distributions of BA corresponding to attractors whose amplitudes change/remain constant after the response system is driven are shown in Figures12(a) and 12(b), respectively. According to the type and amplitude variation of the attractor, BA can be divided into four regions marked with different colors. Partitioning of BA is shown in Table 1. From Figure 12 and Table 1, it can be seen that the amplitude of the base decreases significantly (corresponding to region A) or remains constant (corresponding to regions C and D) under most of the initial conditions and is amplified (corresponding to region B) in only a very few cases.Figure 12
Distribution of BA corresponding to attractors whose amplitudes (a) change or (b) remain constant after the response system is driven.
(a)(b)Table 1
Partitioning of BA.
PartitioningAttractor before being drivenAttractor after being drivenMarkerTypeAmplitudeTypeAmplitudeRegion APeriodicLargeChaoticSmallBlackRegion BPeriodicSmallChaoticLargeRedRegion CPeriodicSmallChaoticSmallGreenRegion DPeriodicLargeChaoticLargeBlueSelect the initial conditions (−1.7, 0.1), (9.3, 8.5), (1.5, −4.5), and (−4.5, 1.5) from different partitions to study the attenuation of line spectra, where PSD is measured in decibels and 1 is used as a reference. As shown in Figure13, when the initial condition (−1.7, 0.1) from region A is selected, the base operates in the large-amplitude periodic attractor when the system is not driven, and the intensity of the characteristic line spectra is 32.65 dB. When the system is driven, the base operates in the small-amplitude chaotic attractor, and the intensity of the characteristic line spectra is 10.9 dB, which is 21.75 dB lower than before. When the initial condition (9.3, 8.5) from region B is selected, after the control is applied, the base migrates from the small-amplitude periodic attractor to the large-amplitude chaotic attractor, and the intensity of the characteristic line spectra increases from 7.309 dB to 33.19 dB. When the initial conditions from regions C and D are selected, the controlled base remains in the chaotic attractor with small amplitude and large amplitude, and the intensity of the characteristic line spectra is 10.99 dB and 33.18 dB, respectively. The favorable initial conditions for keeping small peaks of the line spectra are regions A and C. As shown in Figure 14, the favorable initial conditions for keeping the device's small response are also in areas A and C. In summary, BA can be divided into different regions with significantly different attenuation effects of line spectra and vibration, so as to guide the selection of favorable initial conditions in engineering practice.Figure 13
Phase trajectories and PSD of the base before/after the response system is driven when the initial conditions are from different regions: (a) region A phase trajectories; (b) region A PSD; (c) region B phase trajectories; (d) region B PSD; (e) region C phase trajectories; (f) region C PSD; (g) region D phase trajectories; (h) region D PSD.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 14
Time histories of the device before/after the response system is driven when the initial conditions are from (a) region A; (b) region B; (c) region C; (d) region D.
(a)(b)(c)(d)


## 4.1. Process of Global Analysis

In this paper, a new global analysis method is proposed to study the distribution of BA: Tracking the path of the attractor based on the Poincaré mapping method and the pseudo-arc-length continuation method and reconstructing the phase trajectory of different attractors through the point on the path. In this way, the types of coexisting attractors and their distribution in the phase space can be visualized, providing a priori information for the determination of BA. Select the analysis plane and divide uniformly small grids, and determine the final behavior of the center points of all grids based on the information of the coexisting attractors. The method is described in detail below.In general, anN-dimensional continuous-time nonautonomous system can be transformed into an autonomous one by taking time t as a new one-dimensional variable(10)x˙=fx,s,tt˙=1x,s,t∈RN×RL×R+,where xt∈RN is the state variable, s∈RL is the parameter variable, t is the time variable, and f:RN×RL×R+⟶RN is the vector field. If fx,s,t+T=fx,s,t is satisfied, then a global Poincaré section Σ=x,s,t∈RN−1×RL×R+|t=t0 and Poincaré mapping Px0=xx0,t0+T can be defined. Setting arbitrary initial condition x0,t0 that all phase trajectories of the system intersect with the section Σ transversely, the Poincaré mapping has the same time interval τ=T for all points x0∈Σ. If the Poincaré mapping is continuously differentiable, when the variable parameter changes continuously in small steps, the attractor will move smoothly on the section Σ , so that the path of the attractor can be tracked [24].In this paper, the pseudo-arc-length continuation algorithm is used to track the path of the attractor, which is a prediction-correction method with the following basic ideas: First, determine the variable parameters and iterative step length, and give the preselected parameters and the direction of parameter change. Set the initial conditions; different initial conditions correspond to different Poincaré sections. Then start the prediction-correction process, use the current parameters to predict the steady-state trajectory of the system, and use a convergence algorithm such as the Newton-Raphson method to refine the prediction. The path tracking of the attractor is achieved through an iterative process of modifying parameter values. The path of the attractor under different initial conditions is obtained by changing the direction of the parameter change. It is worth noting that path tracking can only qualitatively analyze the attractor types under different initial conditions but cannot determine the amplitude information. This is because the amplitude of the path is closely related to the selected Poincaré section, so the true amplitude of the attractor cannot be displayed. To this end, we use points on the path to reconstruct the phase trajectories of different attractors and then obtain the distribution of attractors in the phase space.Schematic diagram of the global analysis method is shown in Figure7. The distribution of the coexisting attractors in the phase space is used as a basis for judging the different attractors. The analysis plane P:y1,…yi,yi+3…yn=0,yi+1,1≤yi+1≤yi+1,2,yi+2,1≤yi+2≤yi+2,2 and the restricted state space E:yi+1,3≤yi+1≤yi+1,4,yi+2,3≤yi+2≤yi+2,4 are selected to distinguish different attractors. The analysis plane is finely and evenly meshed. The center point x0 of the initial grid is selected as the initial condition, and the phase trajectories of the attractors are obtained by numerical integration methods such as fourth-order Runge-Kutta. According to whether the phase trajectory is in the restricted state space E, the attractor corresponding to the point x0 is determined. After that, a new center point is selected on the analysis plane to carry out the above analysis process until the final behaviors of all center points can be determined; then all the attractors and their BA in phase space are obtained.Figure 7
Schematic diagram of the global analysis method.


## 4.2. Coexisting Attractors

The model of NVIS established in this paper is a typical nonlinear system with cubic nonlinear stiffness. Moreover, the nonlinear system has a jumping phenomenon. In some intervals of the bifurcation parameters, the same bifurcation value corresponds to three different values of the response amplitude. Among them, the intermediate value corresponds to an unstable state, and the large and small amplitudes that may appear in the mechanical system depend on the initial conditions. In order to analyze the internal mechanism of the coexisting attractors in this system, the path of the attractor is tracked by bifurcation analysis, and the continuation direction is changed to determine the parameter range of coexisting attractors.Poincaré sectionΣ=x,s,t∈RN−1×RL×R+|t=0 and Poincaré mapping Px0=xx0,T are chosen such that all intersections of the phase trajectory with Σ have the same time interval τ=T=2π/ω. With the excitation amplitude f as the variable parameter and taking a step size of 0.01, path changes of attractors for forward continuation (the value of f varies from 0 to 50) and backward continuation (the value of f varies from 50 to 0) are shown in Figures 8(a) and 8(b), respectively. It can be seen that the limit cycle of the system flow is converted into a fixed point on the Poincaré section by Poincaré mapping, and there exists a jump interval f∈15.45,24.37 for the fixed point under different initial conditions. In this parameter interval, there exist two coexisting period-1 attractors for the response system.Figure 8
Path changes of attractors for (a) forward continuation and (b) backward continuation.
(a)(b)Without loss of generality, the pointf = 20 is chosen to reconstruct the phase trajectories of the attractors. As shown in Figures 9(a) and 9(b), when the response system is not driven, set the initial conditions as (0,−5,0,−5) and (0,−5,0,0), respectively, the phase trajectories are two closed loops with significantly different amplitudes, and the power spectra density (PSD) has obvious characteristics of line spectra. It can be inferred that there is a large-amplitude period-1 attractor and a small-amplitude period-1 attractor coexisting in the system; the number and type of the attractors verify the results in Figure 8. As shown in Figures 10(a) and 10(b), phase trajectories of coexisting attractors when the response system is driven are entangled and spread across a certain phase space, and PSD is also converted from line spectra to continuous spectra, which verifies the effectiveness of the control.Figure 9
(a) Phase trajectories and (b) PSD of coexisting attractors when the response system is not driven.
(a)(b)Figure 10
(a) Phase trajectories and (b) PSD of coexisting attractors when the response system is driven.
(a)(b)


## 4.3. Partitioning of BA

The analysis planey1=0,y2=0,−10≤y3≤10,−10≤y4≤10 is chosen, and the small-amplitude attractor is confined in the state space E:−0.05≤y3≤0.4,−0.4≤y4≤0.4, which is used as the basis for judging different attractors. On the analysis plane, y3,y4∈−10,10 is divided into 100 equidistant parts, respectively, and 104 grids with the same spacing are obtained. The center point of each grid is used as the initial condition, and the corresponding phase trajectory is obtained by numerical integration. If the phase trajectory is in the state space E, the behavior of the point is determined to be a small-amplitude attractor; otherwise it is a large-amplitude attractor.When the response system is not driven/driven, distributions of BA are shown in Figures11(a) and 11(b), respectively, where the blue and red regions correspond to the large-amplitude attractor and small-amplitude attractor, respectively. It can be seen that BA corresponding to different attractors have changed, such as the points (-1.7, 0.1) and (9.3, 8.5) marked in the figure.Figure 11
Distribution of BA when the response system is (a) not driven or (b) driven.
(a)(b)The distributions of BA corresponding to attractors whose amplitudes change/remain constant after the response system is driven are shown in Figures12(a) and 12(b), respectively. According to the type and amplitude variation of the attractor, BA can be divided into four regions marked with different colors. Partitioning of BA is shown in Table 1. From Figure 12 and Table 1, it can be seen that the amplitude of the base decreases significantly (corresponding to region A) or remains constant (corresponding to regions C and D) under most of the initial conditions and is amplified (corresponding to region B) in only a very few cases.Figure 12
Distribution of BA corresponding to attractors whose amplitudes (a) change or (b) remain constant after the response system is driven.
(a)(b)Table 1
Partitioning of BA.
PartitioningAttractor before being drivenAttractor after being drivenMarkerTypeAmplitudeTypeAmplitudeRegion APeriodicLargeChaoticSmallBlackRegion BPeriodicSmallChaoticLargeRedRegion CPeriodicSmallChaoticSmallGreenRegion DPeriodicLargeChaoticLargeBlueSelect the initial conditions (−1.7, 0.1), (9.3, 8.5), (1.5, −4.5), and (−4.5, 1.5) from different partitions to study the attenuation of line spectra, where PSD is measured in decibels and 1 is used as a reference. As shown in Figure13, when the initial condition (−1.7, 0.1) from region A is selected, the base operates in the large-amplitude periodic attractor when the system is not driven, and the intensity of the characteristic line spectra is 32.65 dB. When the system is driven, the base operates in the small-amplitude chaotic attractor, and the intensity of the characteristic line spectra is 10.9 dB, which is 21.75 dB lower than before. When the initial condition (9.3, 8.5) from region B is selected, after the control is applied, the base migrates from the small-amplitude periodic attractor to the large-amplitude chaotic attractor, and the intensity of the characteristic line spectra increases from 7.309 dB to 33.19 dB. When the initial conditions from regions C and D are selected, the controlled base remains in the chaotic attractor with small amplitude and large amplitude, and the intensity of the characteristic line spectra is 10.99 dB and 33.18 dB, respectively. The favorable initial conditions for keeping small peaks of the line spectra are regions A and C. As shown in Figure 14, the favorable initial conditions for keeping the device's small response are also in areas A and C. In summary, BA can be divided into different regions with significantly different attenuation effects of line spectra and vibration, so as to guide the selection of favorable initial conditions in engineering practice.Figure 13
Phase trajectories and PSD of the base before/after the response system is driven when the initial conditions are from different regions: (a) region A phase trajectories; (b) region A PSD; (c) region B phase trajectories; (d) region B PSD; (e) region C phase trajectories; (f) region C PSD; (g) region D phase trajectories; (h) region D PSD.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 14
Time histories of the device before/after the response system is driven when the initial conditions are from (a) region A; (b) region B; (c) region C; (d) region D.
(a)(b)(c)(d)


## 5. Conclusion

The chaotification method can be used to weaken the line spectra component of the ship's radiated noise. Aiming at the main obstacle that plagues the engineering realization of this innovative concept, the main work of this paper and the conclusions drawn are summarized as follows:(1)
The dimensionless dynamic equation of the NVIS on a flexible base is derived. A control method similar to generalized synchronization is used to realize the persistent chaotic motion of the system. The maximum Lyapunov exponent, the conditional Lyapunov exponent, and disturbances are introduced to verify the effectiveness of the control method, as well as its stability and robustness in a large range of parameters.(2)
A global analysis method suitable for high-dimensional systems is proposed, which greatly improves the computational efficiency by reducing cyclic operations, and can visualize the types of coexisting attractors and their distribution range in phase space.(3)
The BA is partitioned, and the initial conditions from different partitions are selected to study the attenuation effect of line spectra and vibration, and the favorable initial conditions for realizing the small-amplitude chaos of the system are given. This method may be extended to the global analysis of other controlled systems.



---

*Source: 7149790-2022-02-04.xml*