Generation of equations for shell models of turbulence in the Maple system

The paper describes a technology for the automated compilation of equations for shell models of turbulence in the computer algebra system Maple. A general form of equations for the coefficients of nonlinear interactions is given, which will ensure that the required combination of quadratic invariants and power-law solutions is fulfilled in the model. Described the codes for the Maple system allowing to generate and solve systems of equations for the coefficients. The proposed technology allows you to quickly and accurately generate classes of shell models with the desired properties.


Introduction
The main idea of shell models of turbulence is to construct a chain of ordinary differential equations describing the energy distribution over the scale spectrum in developed turbulence.
The scheme for constructing cascade models is as follows. The axis of wavenumbers is divided into exponentially expanding zones (shells). A sequence of wave numbers arises k n = q n k 0 , where for a suitable choice of the spatial scale k 0 = 1. The denominator of the progression q can be chosen arbitrarily large enough, most often in computational practice q = 2 is taken [1][2][3].
The dynamics of each field X of a turbulent system is described by the collective variables X n (t). The absolute values of X n (t) are interpreted as measures of pulsations of the field X in the ranges of scales with wave numbers from [k n ; k n+1 ). For these variables, a system of equations is drawn up that is qualitatively similar to the initial equations of the turbulent system. Since the equations of hydrodynamics and magnetohydrodynamics are quadratically nonlinear, we are talking about a quadratically nonlinear dynamic system. The coefficients at the nonlinear terms of such a system determine in the model the energy transfer along the inertial interval of scales in developed turbulence. Complex or real variables are used when constructing shell models . In this work, we will only consider real ones.
It is known that, in the ideal limit, there are conserved quadratic quantities in the equations of hydrodynamics and magnetohydrodynamics. Therefore, it is necessary that the equations of cascade models also possess quadratic invariants, which are analogs of these quantities.
Actually, the main mathematical problem in deriving the equations of cascade models is to calculate the coefficients of nonlinear terms that would ensure the existence of certain quadratic invariants.
In this paper, we describe the developed technology for the automated construction of cascade models, in which the calculation of the coefficients of nonlinear interactions is implemented using one of the computer algebra systems (CAS) --the Maple system. Note that Maple is used not for studying ready-made shell models, but for the automated derivation of the equations of these models themselves. Previously, the authors of this paper considered the use of CAS for the automated derivation of shell models of MHD turbulence [4,5]. In this work, we will consider various turbulent systems -«pure» hydrodynamics, magnetohydrodynamics, convective system for conducting and non-conducting fluid.
All the details of the mathematical calculations can be found in the work of the authors [6].

General form of the studied models
Let U n , Θ n and B n be collective variables for the fields of velocity, temperature and magnetic induction, respectively. In the most general form, the considered shell models can be written in the following form.
• Hydrodynamics: where R -is the Reynolds number, f n (t) are functions that simulate the supply of energy to a turbulent system from the outside. Usually, this supply is performed only on the main scale of the system (macroscale) n = 0.
• Convection for non-conducting fluid: where Pr -is the Prandtl number.
• Magnetohydrodynamics: where R m -is the magnetic Reynolds number.
• Magnetoconvection: The indices n, i, and j in the equations can take any integer values. Various cascade models differ in the coefficients of nonlinear interactions S ni j , H ni j , S ni j , L ni j , W ni j , which ensure the existence of different quadratic invariants.
In the quadratic terms present in the general forms of shell models (1 -4), the coefficients of these forms are determined by three indices, the first of which corresponds to the number of the equation, that is, the number of the selected scale shell, and the other two -the numbers of the shells interacting with the selected one, i.e. there is some inhomogeneity in the scales space.
However, in reality, there is scale homogeneity for nonlinear energy transfer within the inertial interval. In this case, the interaction of the i-th and j-th shells with the n-th should depend not on the i and j themselves, but only on the distance from them to the n-th shell in space scale. In addition, all nonlinear terms of the exact equations of hydrodynamics, magnetohydrodynamics, and convection contain the nabla operator. An analogue of this operator in the n-th scale shell is the wave number k n = q n . Therefore, nonlinear terms for the general form (1-4) of cascade models are usually written in a different way: The left-hand of the notation turns out to be more convenient for obtaining equations for the coefficients that ensure the fulfillment of the conservation laws, so we took it as the initial one. The right-hand version of the record corresponds to the concept of scale homogeneity within the inertial interval and qualitatively reflects the differential structure of the nonlinear terms of equations. Note also that the matrices of the coefficients s i j and l i j are naturally assumed to be symmetric.

Constraints on the coefficients of nonlinear interactions
In the ideal limit, i.e. when R = Pr = R m = ∞, and in the absence of external forces, the equations of hydrodynamics and magnetohydrodynamics satisfy a certain set of conservation laws. These conservation laws differ in the case of «pure» hydrodynamics and in the case of magnetohydrodynamics, and also depend on the dimension of the physical space in which the turbulent system [1,3] is considered. The conserved quantities are some quadratic forms of turbulent fields. Therefore, it is required that the shell models, also possess some quadratic invariants of collective variables. These quadratic forms should be formal analogs of the exact conserved quantities [3,6].
The invariants used in this work are given in the Table 1.
The invariance of each of those given in Table 1 of quadratic forms imposes constraints on the coefficients of nonlinear interactions in the form of homogeneous linear equations.
Formally, these equations, as well as the desired coefficients themselves, are infinitely many, however, restrictions are always introduced on the range of interactions in the space Helicity Kinetic energy Kinetic energy

Magnetohydrodinamics Total energy
Total energy

Cross helicity
Cross helicity

Temperature invariant Energy of temperature pulsations Energy of temperature pulsations
n of scales. This leads to a finite truncation of both the number of nonzero interaction coefficients and the number of equations for them, since the rest of the equations degenerate into identities. The specific form of the equations depends on the conservation laws used. Table 2 presents the equations obtained by us for the coefficients of the matrices of nonlinear interactions, grouped by different types of models, in the case of representation of nonlinear terms by the right-hand sides of formulas (5). The derivation of these equations can be seen in the work of the authors [6].
Let us consider one more restriction that naturally follows from wave interactions and is not associated with quadratic conservation laws. This limitation is fundamental in the construction of non-local in scale models, when formally there are interactions of far from each other scales.
From the representation of the equations of motion of the medium in the space of wave vectors, it follows that those and only those waves participate in the nonlinear interaction, from the wave vectors of which it is possible to make a triangle [1,3]. Then it turns out that for the interaction of shells with numbers n, n + i, n + j, the pair of indices (i, j) must lie in some region, the specific form of which depends on q, but the general the shape is similar to the region F shown in Fig. 1. The figure 1 itself shows the shape of the domain F for q geq2. As q decreases, the region expands. The scheme for calculating such a region of pairs of indices of interacting shells is described, for example, in the works [2,3]. In what follows, we will consider precisely the case q ≥ 2.
To ensure the fulfillment of this constraint, it is necessary to require that all coefficients s i j , l i j , w i j , h i j , whose indices do not fall into the region F, would be zero.
From Table 2 and the requirement that the indices fall into F, it can be seen that the construction of cascade models of all types under consideration is reduced to the compilation and Enstrophy conservation: Magnetohydrodynamics Total energy conservation: Cross helicity conservation: Cross helicity conservation: Magnetic helicity conservation: w i, j + (−1) j w i− j,− j = 0. Symmetries: s i j − s ji = 0, l i j − l ji = 0. 3. Convective system (additional condition for temperature) Energy of temperature pulsations conservation: h i, j + q j h i− j,− j = 0.

Scheme for obtaining equations for coefficients in Maple
Let the range of nonlinear interaction in the space of scales (scale nonlocality) be bounded from above by an integer P ≥ 1. In a real turbulent physical system, such a limitation always exists. This finite scale nonlocality determines the real range of i and j. All coefficients s i j , l i j and w i j , at least one of the indices of which is greater than P in absolute value, will be zero. Since all equations for the interaction coefficients are linear homogeneous, the equations, all the terms of which contain an index whose modulus is greater than P, turn out to be identities. It turns out that the introduction of restrictions on the range of interactions distinguishes finite subsystems in the infinite systems of equations obtained above. It is these subsystems that we will form and solve using Maple tools.
Note also that each equation in Table 2 is defined by a pair of indices (i, j), and this equation contains only coefficients with indices (±i, ± j), (±i, ±(i j)) and (± j, ±(i j)). There are no other combinations. Therefore, only those equations for which the indices (i, j) will not necessarily degenerate into identities, satisfy the condition: They lie in the star-shaped region D shown in Fig. 2. Note that in previous works we used a wider area in the shape of a cross [4,5]. Now it has been reduced. In what follows, we need the predicate d(i, j), which selects the region D and is defined by the expression (6). It is also clear that if we choose the pair (i, j) ∈ D, then |i| ≤ 2P, | j| ≤ 2P and |i j| ≤ 2P. Then the sets of the calculated coefficients are naturally represented in Maple, using twodimensional arrays, the indices of which vary from −2P to 2P. Now we introduce the predicate f (i, j), which selects the region F of pairs of indices of interacting shells (Fig. 1). It's clear that Then all the elements of the arrays, the indices of which satisfy the condition: should be zeros. The rest of the elements must be assigned symbolic values, which will then be the components of the vector of unknowns in the system of equations being compiled.
Next, it is necessary to compose expressions for the left side of the equation from Table 2 for the three-dimensional MHD flow for each pair of indices from the region D. In this case, it is advisable to control whether this expression turns out to be zero, since then the equation degenerates into identity. The corresponding Maple code looks like this: Application of the described scheme and code in the considered case of a threedimensional MHD flow for P = 2 gave the following results (we present only nonzero coefficients): The result is a 2-parameter class of cascade models for any q ≥ 2. Of course, for P = 2 it is possible to obtain expressions for these coefficients manually, without making general equations from Table 2, but simply by reducing terms into derivatives of quadratic forms. However, the advantage of the proposed scheme is that it works for any values of the scale nonlocality P, and allows avoiding errors in algebraic calculations.
A well-known property of developed turbulence are power-law spectral laws within the inertial interval. Since collective variables are interpreted as measures of pulsations of the corresponding scales, the manifestation of the spectral laws in shell models will be stationary power-law solutions.
The system of equations for the model coefficients formed and solved in Maple can be supplemented to ensure the existence of stationary power-law solutions. Let us show how this can be done using the example of the model considered above for a three-dimensional MHD flow. Let in the model (3) R = R m = ∞ and f n = 0. We substitute stationary solutions of the form u n = u 0 k α n = u 0 q nα and B n = B 0 k α n = B 0 q nα . Since by rescaling the variables of velocity and magnetic field it is possible to achieve that u 0 = B 0 , we obtain two homogeneous linear equations for the coefficients of the model: It is these equations that must be added to the system for the interaction coefficients. In the above Maple code, replace line 62 with the following snippet: We can add several groups of equations of the (10) type with different α, if we want to ensure the existence of several stationary power-law solutions. It is clear that the addition of such equations can reduce, but not necessarily reduce, the number of independent variables in the overall solution for the interaction coefficients. In particular, a situation is possible when their addition will lead to a system with only zero solutions. This suggests that there is no cascade model with such a set of invariants and spectral laws.
The pulsation energies of the velocity and magnetic field in the n-th scale shell are determined, respectively, by the expressions u 2 n and B 2 n , which for stationary power-law solutions will obviously be proportional to q 2αn . If, however, the spectral law E(k) ∼ k λ holds for the pulsation energy E(k) in a turbulent system, then the energy in the entire n-th scale shell E n ∼ q (λ+1)n [3]. Therefore, the power-law solution with exponent α = (λ + 1)/2 will correspond to the spectral law for energy with exponent λ in the cascade model.
Then, for example, the Kolmogorov spectrum for energy with λ = −5/3 will correspond to α = −1/3. It turned out that adding equations (10) with α = −1/3 does not change the eqref solve solution. Therefore, we can say that the formulas eqref solve define a class of cascade models in which there is a stationary solution corresponding to the Kolmogorov spectrum.
Of course, having ensured the existence of stationary power-law solutions, we cannot guarantee their stability. However, it is no longer possible to elucidate stability in a general formulaic form. It is necessary to fix the numerical values of the parameters and carry out computational experiments with a specific model. Moreover, it is well known that, in shell models, power-law solutions are often unstable and the spectrum is determined by averaging over time Such calculations are no longer related to the subject of this work. However, in principle, the proposed method allows generating models with power-law solutions.

Conclusions
To study turbulence using shell models, software tools are very useful, which would allow you to quickly generate classes of models that satisfy one or another set of conservation laws. The choice of a specific model within these classes can then be specified using additional physical considerations, for example, the existence of a given probability distribution for the interaction of certain shells.
The article describes the developed technique for constructing cascade models, in which the compilation of a system of equations for the coefficients of nonlinear interactions in the model and its exact solution is implemented using Maple systems The technique allows one to vary the dimensions of nonlocality of nonlinear interaction in the space of scales and expressions for cascade analogs of conservation laws and spectral laws. This technique and the programs that implement it solve the problems of fast and error-free generation of entire classes of cascade models.