The equation of motion for the j-th mass (j = 1, 2, … , 2n + 1. For the center mass of the cell in position i, j = 2i.) from the model in Figure 1 is given by m j u ̈ j t + k j − 1 + k j u j t − k j u j + 1 t − k j − 1 u j − 1 t = f j , (1) where u _j (t) is the displacement of the j-th mass m _j , k _j is the constant of the j-th spring, and f _j is the external force loaded on the j-th mass. According to the periodicity of structure, the value of parameters are m j = m a 2 , j = 1 , 2 n + 1 m b , j is even m a , otherwise , (2) k j = 0 , j = 0 , 2 n + 1 k b , j is even k a , otherwise . (3) And so on, the equation of motion for each mass in the structure can be obtained. Under the assumption of harmonic solution, u ( t ) = u ̃ e i ω t , f ( t ) = f ̃ e i ω t , where u ̃ and f ̃ are the amplitudes of displacement and force, and the equations of motion can be written in the matrix form: − ω 2 M U ̃ + K U ̃ = F ̃ , (4) where U ̃ = [ u ̃ 1 , u ̃ 2 , … , u ̃ 2 n + 1 ] T , F ̃ = [ f ̃ 1 , f ̃ 2 , … , f ̃ 2 n + 1 ] T , M = m a 2 m b m a m b ⋱ m b m a 2 2 n + 1 , (5) K = k a − k a − k a k a + k b − k b − k b ⋱ ⋱ ⋱ k a + k b − k b − k b k b 2 n + 1 . (6) For free vibration, the external force is 0 , and Eq. 4 becomes K − ω 2 M U ̃ = 0 (7) The only non-trivial solutions of Eq. 7 correspond to the values of ω ² that satisfy the characteristic equation K − ω 2 M = 0 (8) Then, the natural frequencies and modal shapes can be obtained. Solving Eq. 4 directly when harmonic excitation is loaded gives the amplitudes of displacement response, U ̃ = K − ω 2 M − 1 F ̃ (9)

Individual mass and spring in the lumped-mass model are analyzed to get the transfer matrix. From the relationship of forces and displacement at both ends, and the harmonic vibration assumption, equations of motion are obtained. For mass, u R F R = 1 0 − m ω 2 1 u L F L , (10) or v R = T m v L , (11) for spring, u R F R = 1 1 k 0 1 u L F L , (12) or v R = T k v L , (13) where u _L, u _R, F _L, and F _R are the displacement and the forces acting on the left and right ends, m and k are the parameters of mass and spring, v _L and v _R are the state vectors of left and right ends, T _m and T _k are the transfer matrix of mass and spring, respectively.

The unit cell in Figure 1 is composed of masses and springs in series; thus, the transfer matrix of unit cell T _uc is the product of transfer matrix of each component expressed as v R = T uc v L , (14) T uc = T m a 2 T k b T m b T k a T m a 2 . (15) According to the Bloch theorem, when a free wave travels in the periodic structure, the following condition should be satisfied: v R = λ v L (16) where λ = e ^ikd is a Floquet multiplier, k is wavenumber, and d is lattice constant. Substituting Eq. 16 in Eq. 14 leads to the following eigenvalue problem, T uc − λ I v L = 0 (17) By solving this eigenvalue problem, the eigenvalues appear in pairs of ( λ , 1 λ ) representing a positive- and a negative-going wave along the propagating direction and the eigenvectors correspond to wave shapes. The eigenvalues with ‖λ‖ = 1 and ‖λ‖ ≠ 1 correspond to propagating waves and evanescent waves, respectively (Fan, 2016). Actually, for the case of the 1D diatomic lumped-mass model, the expressions of the boundaries of the band gap can be obtained easier by using the equations of motion for two neighboring masses and the Bloch theorem directly (Hussein et al., 2014). Also, the range of the first Bragg band gap is m a + m b k a + k b − Δ 2 m a m b , m a + m b k a + k b + Δ 2 m a m b , (18) where Δ = m a + m b k a + k b 2 − 16 m a m b k a k b . (19)

For 10⁴ samples with α = 10%, the natural frequencies are calculated by Eq. 8 and shown in Figure 3. The nominal natural frequencies are marked by diamond, and the bars are the range of variation of natural frequencies for the samples. The 1st Bragg band gap of the model without disorder is also obtained as a reference by Eq. 18, which is (4.47 and 6.32) consistent with the result from Figure 2. It can be seen from Figure 3 that when disorder does not exist, there is no natural frequency in the Bragg band gap, and the 11th natural frequency is coincident with the stop frequency of the band gap. Because this finite structure with free ends is symmetric, the natural frequencies are within or at the bounds of each propagation zone (Mead, 1975). Moreover, the influence of disorder on natural frequencies at the beginning of each propagation zone can be almost ignored but becomes gradually important with the increase of order. The range of variation of the 10th natural frequency, which is close to the start frequency of band gaps, has overlapping areas with the Bragg band gap.

Next, we consider forced response and calculate amplitudes of displacement at some specific frequencies within the band gap. The responses are shown in Figure 4, where the red diamond markers are the amplitudes of the nominal displacement response. It can be seen that the amplitude of masses in the same cell decreases with the increase of the number of masses for all the frequencies. This is the typical characteristic of band gaps. The propagation of waves is prohibited within the band gaps exponentially; thus, the energy is confined near the origin of excitation and the structural response tends to be localized to the immediate vicinity of the excitation source (Langley et al., 1997).

We can also notice that the position of maximum amplitude is the second mass for f = 4.7 and 4.9. For the rest of the frequencies it is the first mass. This can be explained from the perspective of wave mode. We calculated the wave mode shapes in the 1st Bragg band gap. On a special note, we call that Mass 1, 2, and 3 are components of the unit cell, not the structure, and Mass 2 is the center mass. When cells are assembled to a periodic structure, Mass 2 corresponds to masses in the structure with an even index, while Mass 1 and 3 do not correspond directly to masses in the structure except the first and last masses. Mass 3 must be with Mass 1 in the next unit cell together to constitute a whole mass. Therefore, the amplitude of Mass 3 is equal to that of Mass 1 in the next cell, of which the amplitude is decided by Mass 1 in this cell according to Eq. 16. So, Mass 3’s amplitude is decided by Mass 1’s in the same cell, and only Mass 1 (or Mass 3) and Mass 2 are independent.

The wave mode shapes at frequencies f = 4.7 and f = 6.1 are shown in Figure 5A, B. These two wave mode shapes change with frequencies. Then, the proportion of magnitudes of wave mode shapes in the band gap is calculated and presented in Figure 5C. It can be seen that the proportion of Mass 2 is prominent at the beginning, while the proportion of Mass 1 is almost zero. With the increase of frequency, the proportion of Mass 2 decreases gradually until to zero, and Mass 1 begins to dominate. Actually, we can find the same phenomenon from the result of forced response in Figure 4. The amplitudes of masses with an even number, which corresponds to Mass 2, are the largest in each unit cell for f = 4.7, 4.9 and turn to be smaller gradually than those of masses corresponding to Mass 1 for f = 5.1, 5.4, 5.7, 6.

For the samples with the disorder, the blue bars in Figure 4 are the ranges of variation of amplitudes. It is obvious that amplitudes of displacement are affected by random disorder, and the distributions of ranges differ in positions and frequencies. It is especially noted that there are always cases, where the maximum upper cap of bars exceeds the maximum red diamond, which means the maximum amplitude is amplified. Compared with the nominal situation, the difference for the maximum amplitude position is that it changes to the first mass for f = 4.9 in rare cases, but still in the first unit cell. In vibration control, the amplification of vibration is expected to be avoided. Hence, it is necessary to describe the variation of maximum amplitude quantitatively and figure out the effects of random disorder on it.

To quantify the influence of disorder on the maximum response in a certain band gap (f _a , f _b ), the frequency-averaged spatially maximum response β ̄ is proposed as follows, β ̄ = 1 f b − f a ∫ f a f b max U d max U o d ω (20)

where ω is the excitation frequency and U _d and U _o are the amplitudes of displacement when the structure is disordered and ordered. This indicator is normalized by the nominal maximum response and averaged in the frequency domain. So, it estimates the overall deviation of dynamic performance in the target band.

In this work, we discretize the frequency band (f _a , f _b ) uniformly by N frequency points and calculate the maximum response by FRFs at each frequency. The frequency-averaged spatially maximum response β ̄ turns out to be β ̄ = 1 N ∑ i = 1 N max U d max U o ω i (21) It needs to be mentioned that the target frequency band does not always coincide with the band gap (f _a , f _b ). The magnitude of the imaginary part of the propagation constant kd reflects the degree of attenuation. At the beginning and ending parts of the band gap, the magnitude is small (Hussein et al., 2014), and so the vibration attenuation is limited. Therefore, it is preferred that the target band is a proper subset of band gaps. In this work, we suppose the target band is (f _a + Δ, f _b − Δ), where Δ = 0.1 (f _b − f _a ).

We consider that disorder strength α varies from 1% to 10%. For each situation, the frequency-averaged spatially maximum response of 2 × 10⁴ samples is calculated. The estimated probability density functions (PDFs) are shown in Figure 6. β ̄ = 1 means the same vibration performance with the nominal structure. β ̄ > 1 means the response amplification and β ̄ < 1 is the opposite. With the increase of α, some important results are noticed. First, the majority values of β ̄ is centered at the nominal value 1. Specially, the means shown in Figure 7A do not change significantly but still increase slowly. Second, the areas between PDFs and the x axis at the two sides of the reference line are almost the same. This is checked by the ratio of areas shown in Figure 7B. So, the probability of β ̄ > 1 and β ̄ < 1 is almost equal. The PDFs become flatter, which is consistent with the increase in standard deviation shown in Figure 7A. Thus, β ̄ greater than 1 originally will be more likely to be bigger. In addition, the range of β ̄ is enlarged. The 99.9th percentile shown in Figure 7A increases even with the increase of disorder strength. Specifically, when α = 10%, the 99.9th percentile nearly reaches 1.1.

Once a periodic structure is manufactured, the disorder must exist, and the stochastic results tell that we cannot foresee the performance variation of this structure. The amelioration and deterioration of performance are both possible. It might perform better or worse with the same probability. But with greater disorder strength, the deviation of vibration performance from the nominal structure is bigger and the worst cases get worse in terms of probability. So, the disorder strength needs to be controlled to a minimum.

In this model, there are 10 independent variables (center mass parameters of 10 cells) distinguished by the position number. To investigate the impact of these variables on the indicator defined before, GSA is used. GSA is a better choice than local sensitivity estimated by the partial derivative. Derivatives are only informative at the base point, where they are computed (Saltelli et al., 2008). However, we need to estimate the importance of variables in an n-dimensional box, not only at a specific reference point. Moreover, the local sensitivity is not able to capture the interaction among parameters (Christopher Frey and Patil, 2002).

Variance-based methods have assessed themselves as versatile and effective for sensitivity analysis and look at the entire factor distribution, using customary Monte Carlo methods of various sophistication modes (Saltelli et al., 2010). Among them, Sobol’ method (Sobol, 2001) is extensively applied to calculate the global sensitivity index. Supposing the model f(x) defined in I ⁿ , where I is the unit interval [0, 1], the variance-based framework is described by the ANOVA-representation, f x = f 0 + ∑ s = 1 n ∑ i 1 < ⋯ < i s n f i 1 ⋯ i s x i 1 , … , x i s , (22) where 1 ≤ i ₁ < ⋯ < i _s ≤ n, and ∫ 0 1 f i 1 ⋯ i s x i 1 , … , x i s d x k = 0 f o r k = i 1 , … , i s . (23) It can be derived from Eq. 23 that the members in Eq. 22 are orthogonal and can be expressed as integrals of f(x). Therefore, by squaring Eq. 22 and integrating, we get D = ∑ s = 1 n ∑ i 1 < ⋯ < i s n D i 1 ⋯ i s , (24) where D = ∫ f 2 d x − f 0 2 , D i 1 ⋯ i s = ∫ f i 1 ⋯ i s 2 d x i 1 ⋯ d x i s . (25) D and D i 1 ⋯ i s are the variances of f(x) and f i 1 ⋯ i s . Equation 24 means the variance of output is decomposed into fractions attributed to inputs. Then, sensitivity index is defined as S i 1 ⋯ i s = D i 1 ⋯ i s D , (26) where the integer s is the order of the index. This index reflects the proportion of variance, which means the importance of impact on the output. The bigger the index, the greater the change of output produced by the variation of corresponding inputs and vice versa.

Under this framework, the ranking of input variables can be determined. In practice, the first-order indices also termed the main effects (ME), and the total sensitivity index (TSI) are computed, thus giving good information on the model sensitivities (Iooss and Lemaître, 2015). For the variable x _j , the ME is S _j . It represents the share of the output variance that is induced by the variable x _j alone. A high ME means that a variable is important, but a low ME does not mean that the variable has no influence, as it can be involved in interactions (Wu et al., 2019). The TSI, S T j = ∑ s = 1 n ∑ i 1 , … , i s ∩ j ≠ ∅ n S i 1 ⋯ i s . (27) contains all the sensitivity indices that x _j participates in, including the first-order index S _j . It is a measure of the share of the variance that is removed from the total variance when the considered variable is fixed to its reference value (Christen et al., 2016). Therefore, variables with low TSI can be considered not influential.

We employ Sobol’ method to calculate the ME and TSI of the 10 inputs m _i (j = 2i, i = 1, 2, … , 10. i is the position number, specifically, Position 1 is the closest to the excitation source, and Position 10 is the farthest.) to the output β ̄ . Sobol’ sequence is used to get 1.5 × 10⁵ samples uniformly distributed over the 10-dimensional box. Sampling based on Sobol’ sequences performs better than other common sampling techniques for most of the analyzed aspects (Burhenne et al., 2011).

Figure 8 shows the results of ME and TSI. There are 10 inputs, which correspond to 10 cells, respectively. It is obvious that the ME and TSI of the input at Position 1 are prominent. Their values all exceed 0.9, while the rest is less than 0.1, and almost 0, when the position number is greater than 2. This means that the model output (frequency-averaged spatially maximum response) is most sensitive to the variation in Position 1. For the input at Position 1, the value of TSI is slightly greater than that of ME, so the main contribution of TSI is the first-order index (ME) and other order indices are little. This means the interaction effect of the input at Position 1 with others is weak. That is to say, if we want to achieve the biggest variation of the output with individual change in a single input, the input at Position 1 will be expected to be the best choice. And supposing we changed the input at Position 1, the variation will not increase, if we keep on changing inputs in other cells. In addition, the ME and TSI decrease sharply with the increase in position number.

These conclusions can be understood in the view of wave propagation. In the band gaps, the maximum energy arises in or near the exciting source, and this energy is little affected by the more remote parts of the system (Langley et al., 1997). Specifically, the amplitude decreases dramatically after the propagation in the first few cells due to the characteristic of exponential decay. That is to say, only the variations in the first few cells matter, and other cells far from the excitation point are not able to influence the maximum energy, which is exactly what the results of GSA reveal. The position, where maximum energy arises is verified and shown in Figure 9. It shows the proportion of positions, where the maximum amplitude occurs in the target band. The first and second masses take up all the proportion and they are both in the cell at Position 1. The difference in proportion between the first and second mass is affected by the wave mode vector illustrated in Section 3.1. From the statistical view, the proportion of samples with the disorder is the same as the nominal situation. It is explained by the maximum position of 200 disorder samples shown in Figure 10. The cut-off frequency of the maximum position varies around that of the nominal case (yellow line). Overall, the maximum amplitude still takes place in Position 1 in the disorder situation.