An intriguing analogy of Kolmogorov’s scaling law in a hierarchical mass–spring–damper model

In Richardson’s cascade description of turbulence, large vortices break up to form smaller ones, transferring the kinetic energy of the ﬂow from large to small scales. This energy is dissipated at the small-est scales due to viscosity. We study energy cascade in a phenomenological model of vortex breakdown. The model is a binary tree of spring-connected masses, with dampers acting on the lowest level. The masses and stiffnesses between levels change according to a power law. The different levels represent different scales, enabling the deﬁnition of “mass wavenumbers.” The eigenvalue distribution of the model exhibits a devil’s staircase self-similarity. The energy spectrum of the model (deﬁned as the energy distribution among the different mass wavenumbers) is derived in the asymptotic limit. A decimation procedure is applied to replace the model with an equivalent chain oscillator. For a signiﬁcant range of the stiffness decay parameter, the energy spectrum is qualitatively similar to the Kolmogorov spectrum of 3D homogeneous, isotropic turbulence and we ﬁnd the stiffness parameter for which the energy spectrum has the well-known − 5 / 3 scaling exponent. Kolmogorov’s


Introduction
Many natural phenomena and engineering processes exhibit energy transfer among a range of different scales. The simplest, frequently studied examples of energy transfer are linear and nonlinear chain oscillators [1][2][3][4][5][6][7] . Dissipation is a key feature of such systems. Jiang et al. [8] and Tripathi et al. [9] investigated and compared the performance of nonlinear energy sink (light mass with viscous damping and an essentially nonlinear spring attached to the primary system) and linear tuned mass damper on multi-degree-of-freedom systems. Felix et al. [10] studied energy pumping and the beat phenomenon in a nonideal oscillator attached to a nonlinear energy sink. Pumhössel found a method to effectively transfer energy among vibration modes using impulsive force excitations [11]. Liu and Jing [12] proposed an X-shaped structure which has advantageous vibrational energy harvesting performance.
Turbulent flow is a prime example of a process which exhibits different scales and an energy cascade. Richardson's so-called eddy hypothesis [13,14] argues that the largest eddies are unstable and break up forming several smaller vortices, gradually transferring the kinetic energy of the flow to smaller scales. The turbulent energy cascade is characterized by the energy spectrumÊ(κ) which describes the distribution of the total energy E = Ê (κ)dκ among the different scales. The wavenumber κ ∼ 1/L is associated with the vortex having characteristic size L. The Kolmogorov spectrum (Fig. 1) illustrates the main features of the turbulent energy cascade. It is the energy spectrum corresponding to 3D homogeneous isotropic turbulence assuming a forcing at the large scales [15][16][17]. The energy input mainly affects the large scales (characterized by L 0 ); hence, the bulk of the energy is contained in the large eddies (energy containing range). In the inertial range (intermediate wavenumbers 1/L 0 < κ < 1/η), the energy spectrum is described by the famous Kolmogorov scaling lawÊ(κ) ∼ κ −5/3 [15,17]. The largest wavenumbers (κ > 1/η) are associated with length scales smaller than the Kolmogorov length scale η. This is the so-called dissipation range in which most kinetic energy is lost due to viscous friction. The Kolmogorov spectrum is confirmed by numerous measurements and simulations [18][19][20][21]. For turbulent flows, the energy flux Π(κ) is also defined which shows the rate of energy transfer at the different wavenumbers κ. For 3D homogeneous isotropic turbulence, the energy flux is constant: Π(κ) = , where = −dE/dt is the dissipation rate of the kinetic energy E. Though energy cascades are usually related to nonlinear phenomena, infinite-dimensional linear systems can also contain a wide range of scales and exhibit com-plex behavior [22]. In fact, nonlinear systems can be embedded in infinite-dimensional linear systems [23].
Motivated by this, we asked the question: Is there a simple linear model whose energy spectrum has features similar to that of the Kolmogorov spectrum? To answer this question in the affirmative, we present a linear phenomenological model of turbulence inspired by Richardson's cascade description. We show that this hierarchical model exhibits the characteristic spectrum of energy transfer seen in turbulent flows. This new point of view will hopefully stimulate research connecting eigenvalue distributions of infinite-dimensional linear systems with energy transfer and resonance capture in nonlinear systems. We also state what this model is NOT: It is not meant to explain details of turbulence; in fact, the behavior of the model is different in many aspects.
Although Fig. 1 corresponds to a process which includes continuous energy injection through forcing, our model freely decays its initial energy. There also exist decaying processes with Kolmogorov-like energy spectrum (see, e.g., [19,24]). The energy decay in our model is exponential as opposed to the power law decrease by decaying turbulence, and the energy transfer is not reversible.
Our phenomenological model of vortex breakdown is a binary tree of masses and springs with dampers at the ground level ( Fig. 2a shows a 3-level model). The masses represent the different scales, the springs are responsible for the energy transfer, and the dampers dissipate energy at the lowest scales. Finite-dimensional modeling of turbulence is common; see, e.g., shell models [17,21] and low-dimensional modeling of turbulence with proper orthogonal decomposition [25]. A similar tree-structured linear mass-spring-damper system was used before to tackle other engineering problems, e.g., viscoelastic behavior [26] and damage detection [27]. It was shown in these papers that the infinite version of their system is described by a fractional (1/2th order) differential equation.
In this paper, we intend to produce an energy spectrum similar to that observed in turbulent flows ( Fig. 1) with a simple linear system. We define the scales, wavenumber, energy spectrum, and energy flux of our phenomenological model analogously to those of turbulent flow. In previous papers, we showed preliminary numerical results for the energy spectrum of decaying and forced systems [28] and an analytical derivation of the energy spectrum of our model [29]. We also per- formed numerical calculations on a nonlinear version of the model [30].
The main finding is that the energy spectrum exhibits the main features of the Kolmogorov spectrum shown in Fig. 1 for suitably chosen parameters. By fine-tuning of the parameters, the midrange of the energy spectrum also has the famous −5/3 scaling exponent. Moreover, we also show how the energy spectrum of the model behaves for different parameters and how this is related to the self-similar structure of the model through the eigenvalues.

Mechanistic turbulence model
Our mechanistic model (Fig. 2a) has n levels, level l ∈ {1, . . . , n} consists of 2 l−1 masses, and the total number of masses is N = 2 n − 1. The top mass and the bottom masses are connected to the fixed ceiling and ground, respectively. Every mass m i has a parent and two (left and right) children whose indices are P(i) = i/2 , L(i) = 2i, R(i) = 2i + 1, respectively (i = 1, . . . , N ), where . denotes the floor operation. The equations of motion are with the boundary conditions (the ceiling and the ground are motionless) where x is the vector of displacements and M, C, and K are the mass, damping, and stiffness matrices, respectively.

Model parameters
We introduce a quantity analogous with the wavenumber of a turbulent scale, the "mass wavenumber" Here M l is the mass scale representing masses in level l (average, for example). For the ease of exposition, we set every mass, stiffness, and damping coefficient equal within a level l and denote these with M l , K l , C l (K l and C l refer to the springs and dampers connecting the masses of the l − 1th and lth levels.). The masses are gradually decreased in lower levels, analogously to the eddy hypothesis, in which a large vortex breaks up into smaller ones. The mass scales M l and the stiffnesses K l are specified by power law distributions, while the dampings are tuned to the stiffnesses (by design, only C n+1 is nonzero), i.e.
Thus, the sum of masses in each level is 1. The characteristics of the energy spectrum of the mechanistic model (defined in Sect. 4) mainly depend on the relation between the base parameters of Eq. (4) (1/2 and σ ). This means that one such base parameter is enough to describe the behavior of the mechanistic model without the loss of generality.

The eigenvalue distribution
The eigenvalues of system (2) are the roots of the characteristic equation P(λ) = det(Mλ 2 + Cλ + K) = 0. Claim The characteristic polynomial of the undamped n-level system (2) (C = 0) with parameter σ = 1/2 has the form (6) where U i denotes the ith Chebyshev polynomial of the second kind [31]. Figure 3 shows the eigenvalue distribution of the purely imaginary eigenvalues of system (2) for the undamped (C = 0) case with 10 levels, σ = 1/2. This eigenvalue distribution exhibits a "devil's staircase"type self-similarity which is related to the self-similar structure of the mechanistic model (binary tree structure). He et al. [32] and Kalmár-Nagy et al. [33] reported different classes of self-similar graphs whose adjacency matrices have devil's staircase-type spectrum. In the latter paper, it was also shown that the characteristic polynomial of the adjacency matrix is a product of Chebyshev polynomials. For damped systems (even for different σ values), the frequency distributions (the imaginary parts of the complex eigenvalues) exhibit similar fractal characteristics as that in Fig. 3.

The energy spectrum
For turbulent flows, the energy spectrum shows the "contribution" of the different scales of eddies to the total energy of the flow, i.e. the mean energy stored in the different wavenumbers κ. An analogous energy spectrum is now defined for the mechanistic model showing the energy fraction stored in the different scales M l (the energy spectrum of the finite-sized mechanistic model is discrete). The total mechanical energy E(t) of the system is The total energy is simply the sum of level energies, i.e.
is defined so that the potential energy of the springs connecting two masses is distributed equally between the two levels and the potential energy of the top and bottom springs is added to E 1 (t) and E n (t), respectively: Here, I (l) = {2 l−1 , . . . , 2 l − 1} are the indices on level l, and δ is the Kronecker delta (δ i, j = 1 for i = j, 0 otherwise). This definition of E l (t) ensures its nonnegativity for t ≥ 0. Figure 4 shows E(t), E 4 (t), E 8 (t) of a 10-level mechanistic model with σ = 1/2 (launched from random initial conditions). The behavior of the energies E(t), E l (t) (l = 1, . . . , n) consists of a transient part and a decaying oscillatory state (the dashed lines in Fig. 4a are exponentially decaying functions with the same exponent). This is the characteristic behavior of E(t) and E l (t) (l = 1, . . . , n) as t → ∞.
The asymptotic behavior of the linear system (2)and so of its energy-is determined by its rightmost pair of eigenvalues λ slow , λ * slow = α ±βi (corresponding to the slowest decaying motion, i.e. rightmost in the complex plane). Figure 4b shows how the real and imaginary parts of λ slow , λ * slow are related to the exponent and oscillation period of E(t), E l (t): The exponent is 2α, while the oscillation period is π/β. We define the mean energyĒ in the asymptotic limit as and we similarly define the mean level energyĒ l for l = 1, . . . , n. The mean energy is simply the sum of the level energies, i.e.Ē = n l=1Ē l . The energy fraction stored in mass wavenumber κ l is defined aŝ and expresses the contribution of a scale M l to the total energy of the system. TheÊ l values constitute the discrete energy spectrumÊ of the mechanistic model.

Derivation of the energy spectrum
A more detailed version of this derivation can be found in [29]. Equation (2) can be recast as a system of 2N first-order differential equationṡ y(t) = Ay(t), y(0) = y 0 with where * and − * denote the conjugate transpose matrix and the inverse of the conjugate transpose matrix, respectively [34]. The solution y(t) of Eq. (11) is The total energy of the system is Since we are interested in the asymptotic behavior of the system, only the terms related to the rightmost eigenvalue pair λ slow , λ * slow = α ± βi of A count in the solution of Eq. (11). All the other terms related to the other eigenvalues of A decay much faster and hence can be neglected. Thus, the general solution of Eq. (11) in the asymptotic limit t → ∞ is where b 1 , b 2 are constant vectors depending on the initial condition y 0 . Substituting Eq. (15) in Eq. (14) leads to where B 1 , B 2 , and B 3 are initial condition dependent constants. The decay 2α and the frequency 2β are due to the quadratic form of (14). To perform the averaging in Eq. (9), we substitute Eq. (16) in Eq. (9). The result is The mean level energiesĒ l have the same form as Eq. (17) with different constants B l,1 , B l,2 , B l,3 (l = 1, . . . , n). The energy spectrum is obtained by substi-tutingĒ andĒ l in Eq. (10), i.e.
l = 1, . . . , n. To calculate the coefficients B 1 , B 2 , B 3 , and B l,1 , B l,2 , B l,3 , we specify the initial condition y 0 in the two-dimensional "slow" subspace of the 2Ndimensional state space of Eq. (11). For example, we can take y 0 = Re s slow , where s slow is the eigenvector corresponding to λ slow . For this initial condition, the evolution described by Eq. (16) is exact, i.e.
for t ≥ 0. The unknown constants B 1 , B 2 , B 3 can be computed from the system of algebraic equations created by taking the first and second derivatives of the RHS of Eqs. (14), (19) at t = 0. We can calculate the initial values E(0),Ė(0),Ë(0) from the initial condition y 0 . The solution of Eq. (20) is The constants B l,1 , B l,2 , B l,3 (l = 1, . . . , n) are calculated similarly. Substituting these solutions in Eq. (18) leads tô Even though we can now compute the energy spectrum, its calculation is only feasible for relatively low number (order 10) of levels due to the rapidly increasing size of matrix A (an additional level almost doubles the number of masses). To overcome this issue, we reduce the mechanistic model by replacing entire levels of masses with single representative masses which leads to an equivalent chain oscillator.

The equivalent chain oscillator
Our approach is similar in spirit to renormalization techniques [35]. Renormalization group theory is extensively used to decimate the elements of large vibrational systems/chain oscillators [36,37]. Figure 2b shows the equivalent chain oscillator which should exhibit the same energy spectrum (the energy distribution among the masses in the chain) as the hierarchical model shown in Fig. 2a. The elements of the mechanistic model are connected in parallel within a level; thus, it can be shown (see Eq. 4) that the parameters of the equivalent chain oscillator are The matrix form of the equations of motion is Figure 5 depicts the eigenvalues of a 10-level mechanistic model (λ i ) and those of its equivalent chain oscillator (ν i ) with σ = 1/2. Claim The characteristic polynomial of the undamped system (24) (C = 0) with n masses and parameter σ = 1/2 is where U n denotes the nth Chebyshev polynomial of the second kind (as in Eq. 6 of the first claim). Figure 5 and our two claims (Eqs. 5, 25) support a claim that ALL eigenvalues of the equivalent chain oscillator are eigenvalues of the corresponding mechanistic model for σ = 1/2, and most importantly, the rightmost pair of eigenvalues is equal (ν slow = λ slow ). Even though we have no formal proof, the eigenvalue and energy equivalence seem to be true for the general, σ = 1/2 case.

The energy flux function
For turbulent flows, the energy flux shows the rate of energy transfer at the different scales of eddies. An analogous quantity is defined for the mechanistic model. Similarly to the energy spectrum, a set of fluxes can be calculated, one for each wavenumber. The energy flux Π l (t) of level l is defined as the rate of energy transfer between levels l and l + 1, and it is calculated as Here Π 0 (t) = 0, since the top mass does not transfer any energy to the ceiling. Furthermore, we define the dissipation rate of the system as It can be shown that the recursive formula (26) leads to Similarly to the mean energy, the mean dissipation rate of the system is also defined in the asymptotic limit: The mean fluxesΠ l are defined similarly for l = 1, . . . , n.
Following the framework of Sect. 4, we scale the mean fluxes with¯ (the mean level energies were scaled withĒ) to yield the scaled energy flux for a given mass wavenumber κ l . TheΠ l values constitute the discrete energy flux functionΠ of the mechanistic model.

Derivation of the energy flux function
We are investigating the asymptotic behavior of the system. To perform the averaging in Eq. (28), we substitute the derivative of Eq. (16) in Eq. (28). The result is The same averaging performed on the level fluxes yields That is, we only need the constants B 1 , B 2 , and B l,1 , B l,2 to calculate the energy flux function. In Sect. 4.1, we showed how these constants are related to the initial values E(0),Ė(0),Ë(0), and substituting Eq. (21) in Eq. (32) leads tô It follows from Eq. (33) thatΠ l >Π l−1 andΠ n = 1.

Results
First we compute and compare the energy spectrum of 10-level mechanistic models and their equivalent  Fig. 6. These figures demonstrate that the energy spectrum of the mechanistic model and that of its equivalent chain oscillator are practically the same regardless of the value of σ . This allows us to use the equivalent chain oscillator to calculate the energy spectrum of even larger systems. The energy spectrum of the mechanistic model is calculated with different σ values using the equivalent chain oscillator formulation for n = 20 levels (Fig. 7).
For σ < 1/2 (Fig. 7a) the energy spectrum is similar to the Kolmogorov spectrum (cf. Fig. 1): The peak is located at a small wavenumber and the spectrum has a negative slope in the intermediate range. Figure 7b shows that the spectrum is almost symmetric for σ = 1/2 (the damping breaks the symmetry), and most energy is concentrated in the intermediate scales. For σ > 1/2 (Fig. 7c), the peak is further shifted toward large wavenumbers and the spectrum has a positive slope in the intermediate range. As σ is further increased, local minima appear in the spectrum  Fig. 7d, e). When σ is set even larger, the spectrum changes once more, and the shape depicted in Fig. 7f becomes typical.
In the σ < 1/2 case, we approximate the intermediate range of the spectrum with the scaling lawÊ ∼ κ γ (Fig. 7a, solid line). The main finding is that parameter σ can be tuned to match the Kolmogorov exponent: At σ ≈ 0.492, the exponent is γ = −5/3, agreeing with the Kolmogorov exponent for 3D homogeneous isotropic turbulence.
It was also investigated how the energy spectrum behaves as the number of levels n changes. In Fig. 8, energy spectra for n = 10, 15, 20, 25 are depicted with different σ 's.
For σ < 1/2 (Fig. 8a), the energy spectrum gets extended to include the higher wavenumbers. For σ ≥ 1/2 (Fig. 8b, c), the change of the spectrum is more spectacular: The shape remains qualitatively the same, but the whole spectrum gets stretched to include the higher wavenumbers. This is not just an extension to higher wavenumbers, since the peak is also shifted significantly toward the larger wavenumbers. In the σ > 1/2 case, it is also possible that the shape of the spectrum changes. For example, in Fig. 8c the n = 25 case exhibits a different spectrum shape compared to the other three cases. This change is similar to how the spectrum changes as σ is increased (see Fig. 7c, d).
The changes in the shape of the energy spectrum for the σ > 1/2 cases (Fig. 7c-f) are due to the change in the order of the rightmost eigenvalues and the accompanying significant change of the eigenvector s slow . The same happens when a σ value is given and n is increased as shown in Fig. 8c. In Fig. 9, the real part of four eigenvalues (denoted by λ 1 , λ 2 , λ 3 , λ 4 ) of the 20-level mechanistic model is plotted as the function of σ . This figure illustrates how these eigenvalues-together with their complex conjugate-"compete" with each other to be the rightmost eigenvalue λ slow . Since rightmost refers to the location of the eigenvalue in the complex plane, the rightmost eigenvalue has the largest real part. For σ ≤ 1/2, the rightmost eigenvalue is λ slow = λ 1 . As . Figure 9b, c shows the boxed parts of Fig. 9a in detail.
One can see that λ 2 , λ 3 , and λ 4 become rightmost one after the other. The real part of the eigenvectors corresponding to these four eigenvalues is depicted in Fig. 10 for σ = 0.503 (the equivalent chain oscillator was used here). These eigenvectors have a very different shape. This is why the energy spectrum changes significantly, when the rightmost eigenvalue changes. We also computed the energy flux function of different mechanistic models. Figure 11 shows energy flux functions for different n's and σ 's. As expected, the energy flux functions are monotonously increasing. The maximum isΠ n = 1 which simply means that the energy flux between the last level and the fixed ground is the dissipation rate of the system. Figure 11a shows the energy flux functions corresponding to the σ < 1/2 case. In this case, the energy spectrum is Kolmogorov- like (see Figs. 6a, 7a, and 8a), while the energy flux function is nearly constant (slightly increasing) from the intermediate range to the largest wavenumber. This is as intriguing as the shape of the energy spectrum itself, since the energy flux is constant in 3D homogeneous isotropic turbulence through the inertial range.

Conclusions and future work
Inspired by Richardson's eddy hypothesis, we introduced a mechanistic model of turbulence with different mass scales. The discrete energy spectrum of the model was defined as the fraction of the total energy stored in the different mass scales. A formula (Eq. 22) was derived to calculate this spectrum in the asymptotic limit. Similarly, the discrete energy flux function of the system was defined as the set of energy fluxes between levels scaled by the dissipation rate. A decimation procedure was applied to yield an equivalent chain oscillator with energy spectrum identical to that of the mechanistic model. This enabled the calculation of the energy spectrum and energy flux for many-level systems. In the σ < 1/2 case, the energy spectrum is qualitatively similar to the Kolmogorov spectrum of 3D turbulence. Moreover, σ can be tuned so that the exponent γ of the energy spectrum precisely matches with that of the Kolmogorov spectrum. Furthermore, in this case the energy flux function is practically constant through the intermediate range.
This work is expected to motivate new studies connecting eigenvalue distributions of infinite-dimensional linear systems with energy transfer and resonance capture in nonlinear systems. In future work, the model will be extended to include nonlinearities (e.g., nonlinear energy sinks at the bottom level instead of linear tuned mass dampers), cascade asymmetry (mistuned masses), and instability at the large scales (e.g., negative damping could be added to the system). We also want to study the system with continuous energy pumping at the large or some intermediate scales (this latter would correspond to the case of 2D isotropic turbulence [38,39]) via external forcing.