Influence of coulomb damping on wave propagation behaviors of nonlinear nonconservative phononic chains/lattices

Fascinating nonlinearity-induced behavior of phononic crystals (PCs) has recently become a hot research topic in the community. However, due to the limitations in the analytical modelling of damping in dynamic systems, the study of damped PCs has not received proper attention. In this paper, the influence of Coulomb damping on the wave propagation behavior of cubically nonlinear monoatomic phononic chains is investigated. To do so, the nonlinear dispersion relation is obtained analytically using the well-established multiple scales method and the band structure of the damped nonlinear chains is compared to the ones corresponding to the linear and nonlinear undamped chains. Due to the coupling between the amplitude and the frequency, stemmed from the nonlinear nature of the chain, Coulomb damping can lead to lower dispersion frequencies in the chain. The formulation and results are then expanded to 2D nonlinear lattices. The present manuscript is the first attempt to capture the effect of Coulomb damping on the wave propagation behavior of nonlinear lattices and the results put us one step closer to developing a comprehensive analytical model for the behavior of damped PCs which can in turn lead to invaluable design concepts for nonlinear nonconservative wave-manipulation devices.


Introduction
Wave propagation characteristics of phononic crystals (PCs) have been the subject of numerous valuable studies in the recent years. Their ability to stop the propagation of waves in certain frequencies, which is also known as 'stop-bands', has brought about numerous applications for PCs, including vibration attenuation [1], sound suppression [2,3], impact mitigation [4], energy harvesting [5] and waveguiding [6]. Changing the geometrical formation [7][8][9], material composition [10] or manufacturing process [11] of such structures, their wave attenuation properties can be modulated. This adds an extra layer of tunability to designing wave-attenuation devices and structures [12]. Such characteristics have opened up new avenues for investigating and proposing new mechanisms for wave propagation controls, over the past years [13,14]. Some of PC applications are illustrated in figure 1. Interested scholars are encouraged to study review articles [15][16][17] for more information regarding the various applications of PCs.
While the efforts of the community in reaching a comprehensive model for predicting the behavior of PCs have reached many peaks in the recent years, a lot of subjects in the field have remained unknown up to now. One of the less explored areas in modelling and investigating PCs is the nonlinear structures and the extraordinary phenomena induced by the nonlinear nature of the problem [18]. It was previously shown that unlike linear PCs, in nonlinear periodic structures, the frequency and amplitude are coupled [19]. This exceptional observation is the source of many subsequent studies with the aim of exploring the effects of nonlinearity on the wave attenuation performance of PCs. Investigations in the field include nonlinear monoatomic [20], diatomic [21,22], triatomic [23], mass-in-mass [24,25] phononic chains, chains with coupled translation-rotation coupling [26] as well as continuous structures [27], which are subject to unique behavior due to their nonlinear nature. It was found that the introduction of weak nonlinearity can lead to the exhibition of a new stop-band mechanism beyond the two well-known mechanisms of Bragg scattering and local resonance [28]. In addition, unique phenomena such as wave-wave interactions [29], nonreciprocity [30][31][32], subharmonic stop-bands [33], etc are some of the direct implications of considering nonlinearity in PCs. These nonlinearity-induced phenomena have led to the proposition of more application-centered concepts, such as nonlinear frequency isolators [34] and acoustic diodes [35]. Exquisite behavior of nonlinear phononic crystals is summarized in figure 2.
Besides the 1D chains, 2D nonlinear lattices are also associated with many intriguing nonlinearity-induced phenomena, such as the existence of solitons [36] and amplitude-dependent directionality [37], which can broaden our vision on the effect of nonlinearity on the dynamic behavior of lattices, as well as providing additional means for controlling wave propagation and directionality in such structures. However, compared to the one-dimensional chains, the 2D lattices have received much less attention. More specifically, 2D structures  with out-of-plane degrees of freedom have shown great potential in real-life experimental efforts. Due to the recent advancements in the micro/nano engineering, such structures have gained profound importance in measurements in small-scale structures [38]. Due to the very high frequencies of vibration in small-scale structures, the video-based measurement devices cannot always capture the in-plane motion of micro/nano sized structures. Therefore, due to such practical shortcomings of the current measurement devices, the out-ofplane motion sensing started popularizing. However, lack of analytical and computational investigations on the out-of-plane wave propagation in lattices and grids is quite obvious in the open literature and more thorough study of their behavior is necessary. While there are a number of studies on the numerical and analytical studies of the wave propagation in micro-lattices with out-of-plane displacements, they are mainly focused on the linear structures. As a result, there is a gap in the investigation of damped planar lattices, which we aim to fill in the present manuscript.
In the previous years, the experimental investigations on the wave propagation of PCs have faced a great leap, both in quantity and quality. However, analytical and computational studies of such structures could not keep up with the extent of experimental progress. One such subject in the analytical study of PCs is the treatment of damping in the models. In the literature, many profound experimental research efforts can be found. However, compared to the experimental studies of PCs with various types of damping, the analytical study of damping in periodic structures seems underdeveloped. In previous years, a number of studies considered damping in PCs and its effect on the wave propagation and attenuation performance. In addition to discrete 1D and 2D massspring models [39], the damping effect on the dispersion curves of continuous structures has also been previously investigated [40]. However, the abovementioned efforts mostly describe the dynamic performance of linear PCs and neglect the effect of nonlinearity. As such, the number of studies investigating the nonlinear damped PCs is far lower than the efforts on their undamped counterparts. However, the intertwined effect of nonlinearity and damping is shown to be of great value in modelling and tuning the wave propagation behavior of PCs. Previous studies have investigated the effect of viscous and quadratic damping on monoatomic chains [41] as well as the effect of fractional damping on the wave propagation in nonlinear 1D and 2D monoatomic lattices [42]. Additionally, the analytical study of 1D diatomic nonlinear damped PC is another major contribution to the field [43].
The dry friction damping mechanism dissipates the energy of a system in the form of heat as a result of the sliding motion between two contacting surfaces. As a result, modelling this type of damping is crucial for describing any structure with contact surfaces. However, due to their nonlinear nature, inclusion of Coulomb damping mechanisms in the governing equations of the motion makes the dynamic analysis of the system more complicated. Therefore, it is often neglected in modelling dynamic systems. However, the Coulomb damping cannot always be neglected. Hence, it is necessary to propose a model for predicting the behavior of structures with Coulomb damping. In the nonlinear dynamics community, however, the analytical modelling of systems with Coulomb has received noticeable attention. Examples include studying the chaotic vibration of a nonlinear oscillator with Coulomb damping [44], perturbation-based analysis of a van der Pol system with Coulomb damping [45] and investigating the simultaneous effect of nonlinearity and Coulomb damping on the response of a nonlinear energy harvester [46], among others.
In this research, the wave propagation of weakly nonlinear one-and two-dimensional PCs with Coulomb damping is studied. The weak nonlinearity is assumed to be of the cubic type and a discontinuous expression for the Coulomb damping force is utilized. The method of multiple scales is used for obtaining the analytical dispersion relations for 1D and 2D PCs and investigating the effect of Coulomb damping on the wave propagation behavior of such structures. The results of the present manuscript can significantly facilitate the modelling and predicting the behavior of real-life acoustic metamaterials and be considered as a big initial step towards designing nonlinear PCs with tunable properties.

Formulation
2.1. One-dimensional nonlinear chain with coulomb damping A schematic view of a one-dimensional phononic chain with Coulomb damping is presented in figure 3(a). The nonlinear PC consists of an infinite number of masses m, which are connected to their neighbors by means of nonlinear springs, whose force-displacement relation is described byẽa = + f ku u 3 [34,43]. k andã are the linear and nonlinear spring constants, respectively, and e is a very small book-keeping parameter, which is a measure of the strength of the nonlinearity [47]. As a results, the springs exhibit a weakly nonlinear behavior. In addition, a Coulomb damper is connected to each mass of the chain and exerts a forcef on the mass. The Coulomb damping force is described by equation (2). Also, the red shade in figure 3(a) designates the nth unit cell of the chain. As observed in the unit cell, the displacement of the nth mass is denoted by u , n whileu n 1 and + u n 1 designate the displacement of its adjacent masses.
For the chain described above, the equation of motion of the nth unit cell is described as: In equation (1),ü n is the acceleration of the nth mass of the chain and f denotes the Coulomb damper force on the nth mass, defined as: in which, u n  denotes the velocity of the nth mass and m is the Coulomb damping coefficient. The damping force is illustrated in figure 3 By introducing the non-dimensionalized time t w = t, n with w = / k m , n the non-dimensional form of the governing equation can be expressed as: In the present paper, the method of multiple scales is utilized to find the analytical dispersion relation of the provided chain. More details regarding the multiple scales methods is provided in appendix A. Using this technique, the solution u n can be expanded as: denoting the fast and slow time scales, respectively. Therefore, the second time derivative is expressed with respect to different time scales as: in which, D 0 and D 1 are derivatives with respect to T 0 and T , 1 respectively. Substituting equations (4) and (5) in the non-dimensional governing equation (equation (3)), the first two perturbation orders of the governing equation of motion are expressed as:  The solution of the zeroth order governing equation (6) can be obtained by admitting the plane wave solution as: Moreover, k denotes the wavenumber. Substituting the solution (8) in equation (6), the linear frequency of the chain is expressed as below: Hence, the ( ) e 1 order solution can be expressed as: The solvability condition of equation (10) requires the coefficients of ( known as the secular terms, to be zero. Otherwise the solution u n would become unbounded as the time increases. Therefore, in order to remove the secular terms, the following equality should be satisfied: By expressing the amplitude in the polar form T 1 i 1 with a and b denoting the real-valued amplitude and phase, respectively, equation (11) can be expressed as: In the previous equation, f r and f i designate the real and imaginary parts of the Coulomb damping term, respectively. Moreover, the damping term f can be described as: For more information regarding the procedure of obtaining equation (13), check out appendix B.
Then, separating the real and imaginary terms, one reaches at: Solving the ordinary differential equation (14), one can find the amplitude relation as: Substituting the amplitude a from equation (16) in equation (15) and assuming the initial condition of b = 0 at = T 0, 1 the nonlinear correction frequency of the chain can be found as: which will lead to the following relation for the nonlinear frequency:

Two-dimensional nonlinear lattice with coulomb damping
The schematic view of the square lattice with Coulomb damping is illustrated in figure 4. Each mass is connected to two neighboring masses in a 1 direction and two in the a 2 directions. The springs in the a 1 direction have the linear and nonlinear stiffness k 1 andã , 1 respectively; while the linear and nonlinear spring coefficients in a 2 direction are denoted by k 2 andã , 2 respectively. Moreover, a Coulomb damper is connected to each mass and provides the energy dissipation for the structure. The system has only one degree of freedom, u , m n , which is perpendicular to the plane of periodicity and in the a 3 direction.
For the mass on the nth unit cell, the governing equation of motion is written as Similar to the 1D case, the non-dimensional form of the governing equation can be expressed as: with the parameters defined asã a = /k , 1 1 1ã a = /k 2 2 1 and = / R k k .

1
The displacement u m n , is expanded in the following form: Substituting the displacement u m n , from equation (21) in the non-dimensional governing equation (20) and separating the ( ) e 0 and ( ) e 1 orders, the following relations can be obtained: The governing zeroth order equation (22) can be solved by assuming the following expression for the displacement u :  Therefore, the linear frequency can be obtained as: Similarly, assuming the amplitude to be of the polar form a T e 1 2 T 1 i 1 and separating the real and imaginary parts, the differential equations for finding the amplitude and frequency are described as: Solving the differential equation in equation (28), the amplitude can be obtained as: Therefore, similar to the 1D case,

Results and discussion
In the current section, the results of the wave propagation analysis in 1D and 2D nonlinear lattices with Coulomb damping are presented and analyzed in more detail.

1D nonlinear chain
The time-and wavenumber-dependent amplitude of a nonlinear PC with Coulomb damping is presented in figure 5(a). In nonlinear PCs, the amplitude not only depends on the time, but it depends also on the wavenumber. According to equation (30), the variation in amplitude is linear in terms of time (T 1 ) and nonlinear in terms of wavenumber (k). This can, in turn, lead to time-dependent dispersion curves for the damped nonlinear PC. Figure 5(a) illustrates the variation of the amplitude in the first Brillouin zone ( k p < < 0 ), as the normalized time is passed from = T 0 1 to p = / T 10 . 1 It is clear that the amplitude decays linearly as time passes. On the other hand, the more complex relation between the amplitude and the wavenumber is clearly visible in figure 5. As the wavenumber progresses inside the first Brillouin zone, the amplitude increases. As a result, the longer wavelengths; i.e., smaller wavenumbers, are associated with lower amplitudes in nonlinear PCs with Coulomb damping. Also, the growth in the amplitude is more abrupt at the beginning of the Brillouin zone. However, it slows down as the wavenumber increases along the Brillouin zone.
In order to elaborate the effect of Coulomb damping on the time-dependent dispersion curves of nonlinear 1D phononic chains, the deviation of the frequency of the damped nonlinear chain from the nonlinear undamped chain is illustrated in figure 5(b) in the Brillouin zone of the chain for the normalized time range of p < < / T 0 10. 1 Note that the nonlinear frequency of the damped chain is described by equation (31), while the undamped nonlinear frequency is expressed as . Interestingly, although the amplitude is lower for the longer wavelengths along the Brillouin zone, the frequency deviation is seen to be lower compared to the shorter wavelengths. Also, the decay in the frequency is also visible as time elapses. It is noted that time has a more prominent effect on the frequencies for higher wavenumbers. For better perception of the notion of time-dependent wave propagation in nonlinear PCs, the band structures of a weakly nonlinear chain at three different normalized times are illustrated in figure 6(a). It should be noted that the normalize time = T 0 1 denoted the nonlinear dispersion curves corresponding to the undamped chain. Figure 6(a) clearly shows that passing through time leads to lower frequencies in structures with hardening nonlinearity and Coulomb damping. As a result, the dispersion curves lie further from the undamped frequency and formed in lower frequencies. Consequently, the undamped frequency serves as the upper bound for the values of frequency, while the frequency of the linear chain indicates the lowest value the frequency can attain. Moreover, figure 6(a) illustrates that as the wavenumber is increased from k = 0 to k p = along the Brillouin zone, the time-dependency of the dispersion curves become more prominent.
To elaborate this matter, the deviation of the damped frequency from the undamped one along the Brillouin zone is illustrated in figure 6(b) for different values of normalized time T . 1 Two sets of conclusion can be drawn from this illustration. First, wavenumber plays a crucial role in determining the time-dependency of the nonlinear frequency of the chain with Coulomb damping. More specifically, while the effect of time on the band structure is negligible at lower wavenumbers and the dispersion curves corresponding to all values of time are formed close to each other, they begin to diverge as the wavenumber is increased. for k p = /2. It can clearly be observed in figure 6(b) that the difference between B and ¢ B is higher than the difference between A and ¢ A . Therefore, the time-dependency of the nonlinear frequency is more prominent in higher wavenumbers. Second, the deviation of the dispersion curves of the nonconservative chain from the nonlinear undamped case ( = T 0 1 ) can also be clearly illustrated in figure 6(b). It should be noted that due to the hardening nature of the nonlinearity, the damped dispersion curves are formed below the nonlinear undamped band structure. Also, to rule out the effect of the damping coefficient, a constant value of m = 0.3 is assumed in figure 6.
Moreover, the effect of Coulomb damping coefficient on the dispersion curves of 1D nonlinear PCs is illustrated in figure 7(a). The solid black curve corresponds to the band structure of a nonlinear PC without damping, i.e., a ¹ 0 and m = 0, while the dotted purple curve denotes the dispersion curves of a linear undamped PC (a = 0 and m = 0). It is observed that increasing the damping coefficient in a nonlinear PC leads to more damping and therefore, lower frequencies. As a result, based on the values of the damping coefficient, the dispersion curves of a damped PC can fluctuate between the ones corresponding to the undamped nonlinear and linear band structures.
To better elaborate the influence of Coulomb damping on the dispersion frequency in a nonlinear chain, the deviation of the nonlinear frequency in a damped chain from the nonlinear undamped frequency is presented in figure 7(b) for different values of damping coefficient. The blue dashed curve denotes the nonlinear undamped chain. While the effect of wavenumber on the deviation is similar to figure 6, structures with more damping are subject to more deviation from the nonlinear frequency. Furthermore, the effect of damping is more noticeable in higher wavenumbers, as the dispersion curves corresponding to different damping coefficients are formed closer at the beginning of the Brillouin zone and tend to diverge in higher wavenumbers.
Furthermore, the amplitude ( ) a T 1 and frequency ( ) b T 1 are also obtained from equations (16) and (17).
Compared to the motion of the undamped chain (m = 0), increasing the damping coefficient leads to lower peak amplitudes in each cycle of the motion. The amplitude reduction is observed to be linear, as obvious from equation (16). The linear decay of the amplitude was previously reported for nonlinear vibrational systems with Coulomb damping. It is also evident from figure 7(c) that increasing the damping coefficient reduces the frequency, as the wavelength is seen to be larger in structures with higher damping coefficients. In order to consider the damping coefficient effect only, a constant wavenumber of k p = 3 was adopted to obtain the timedependent amplitude resulting in the time history of figure 7(c).

2D nonlinear chain
The dispersion behavior and dynamic properties of a nonlinear lattice with Coulomb damping is expectedly more complex than a nonlinear chain, due to the more sophisticated definition of the Brillouin zone in 2D lattices. More specifically, unlike the 1D chains, the Brillouin zone in 2D lattices consists of three separate sections (GX, XM and MG) with different orientations. As a result, one cannot expect the effect of wavenumbers on the amplitude and dispersion curves to be monotonic. This is elaborated in more detail in figure 8, in which, the time-and wavenumber-dependent amplitude in a nonlinear chain with Coulomb damping is presented. It is demonstrated in figure 8 that the amplitude in a 2D lattice with Coulomb damping decays in time linearly. Also, the wavevectors in paths GX and MG along the Brillouin zone is shown to have a more noticeable effect on the amplitude, especially in the regions closest to the long-wavelength limit. Moreover, the effect of time and wavevector on the nonlinear frequency of a lattice with Coulomb damping is presented in figure 9(a), where the deviation from the nonlinear undamped frequency is presented in more detail. Figure 9(a) proves that the variation of the nonlinear frequency along the edges of the Brillouin zone is not monotonic and the magnitude and direction of the wavevectors have a significant effect on the frequency. More specifically, while the deviation of the nonlinear frequency attains its minimum values around the high symmetry point G, the deviation tends to increase as the wavevectors are defined along the XM part of the Brillouin zone.
Also, the time-dependent dispersion curves of a nonlinear 2D lattice with Coulomb damping is presented if figure 9(b) for different values of damping coefficient m. For more elaboration, the dispersion curves corresponding to the linear and nonlinear undamped lattices are also depicted as purple dotted and blue solid curves, respectively. Similar to the 1D case, the two abovementioned conditions serve as the upper and lower bounds for the dispersion curves, i.e., depending on the values of the damping coefficient, the dispersion curves are always formed higher compared to the linear dispersion curves and lower than the nonlinear undamped curve. Furthermore, increasing the value of the damping coefficient will move the dispersion curves further from the undamped nonlinear curves and closer to the linear ones.

Numerical validation
In order to validate the semi-analytical formulations and dispersion curves, a finite chain made of 500 masses is numerically modelled and the governing equation (1) is solved for this chain. In doing so, each individual mass of the chain is subjected to initial displacement ( ) where j is the mass number, and w designates the frequency corresponding to wavenumber k. The simulation is performed for a time t = t 5 n and t p w = / 2 n is the oscillation period. The numerical estimation of the points on the dispersion curve is then performed by using fast Fourier transform (FFT) on the spatial and temporal profiles provided by the numerical solution. However, to rule out the wave distortion at the two ends of the chain, the FFT is performed on a sub-chain in the middle of the primary chain; e.g. masses 200 to 400, which is far away from the edges. Comparison of the analytical and numerical dispersion curves of a 1D nonlinear phononic crystal with Coulomb damping is presented in figure 10 and an excellent agreement is observed between the two results.
The extended numerical validations considering the effect of nonlinearity, damping ratio and the bookkeeping parameter e are presented in appendix C.

Conclusion
In the present manuscript, the wave propagation in nonlinear PCs under Coulomb damping was investigated. The nonlinear dispersion equations were found analytically using the multiple scales method and the effect of Coulomb damping on the dispersion curves and wave attenuation behavior was investigated analytically. It was found that in weakly nonlinear chains with cubic nonlinearities, the amplitude is both time-and wavenumberdependent and decays linearly with time and nonlinearly with wavenumber. The decay in the amplitude owing to the Coulomb damping can lead to the formation of dispersion curves at lower frequencies. This is due to the coupling between the frequency and the amplitude in a nonlinear PC, which leads to extraordinary amplitudedependent behavior in nonlinear PCs. Moreover, the formulations were also expanded to 2D nonlinear lattices with out-of-plane degree of freedom. Similar time-and wavevector-dependent behavior is observed inside the irreducible Brillouin zone of a 2D lattice, with lower dispersion frequencies for damped lattices. The semianalytical results were then validated numerically.
Results of the current paper can expand our vision regarding the nonlinear nonconservative PCs and, in turn, enable us to present ideas regarding their application in real-life devices. It is among the first attempts to analytically investigate damped PCs, which can effectively tackle the problem of modelling damping in PCs.

Data availability statement
All data that support the findings of this study are included within the article.

Appendix. A
Multiple scales method Perturbation techniques are semi-analytical approaches to find the closed form solutions for the responses of a weakly nonlinear system by perturbing the response of the corresponding linear system [48]. The multiple scales method is on of such techniques which can provide a uniform expansion of the solution. The main advantage of the multiple scales method compared to other perturbation techniques such as the Lindstedt-Poincare method is that it can be utilized to predict various nonlinear resonance phenomena and damped systems conveniently [47].
Unlike the Lindstedt-Poincare technique, the multiple scales method considers the expansion representing the response to be a function of multiple independent time scales, instead of just a single one. This is due to the fact that the response (u) depends on a combination of time (t ) and the book-keeping parameter (e) rather than t and e, individually. In other words:ˆ( As a result, the initial ordinary-differential equation describing the motion is replaced by a partial-differential equation in terms of different time scales. More information regarding the perturbation methods and the multiple scales method can be found in [47,48].

Appendix. B
The Coulomb damping term f can be expanded in the form of Fourier's series by adopting the new variable ( ) In this way, the Coulomb damping force is described by: