Quentin Pikeroen

Amaury Barral

Guillaume Costa

Ciro Campolina

Alexei Mailybaev

Berengere Dubrulle
email: berengere.dubrulle@cea.fr

Tracking complex singularities of fluids on log-lattices

come

Introduction

Viscous fluids dissipate mechanical energy into heat due to the first law of thermodynamics. Observations and numerical simulations reveal that this dissipation is not homogeneous within the flow but occurs via spatially or temporally intermittent bursts, a phenomenon classically referred to as intermittency. Moreover, after spatial and temporal averaging, the mean energy dissipation becomes independent of the viscosity in the inviscid limit, according to the empirical "zeroth law of turbulence". Onsager explained these observations in 1949 [START_REF] Onsager | Statistical hydrodynamics[END_REF], conjecturing that strong enough singularities in the inviscid flow could provide an anomalous dissipation. While this conjecture has been proven mathematically [START_REF] Isett | A proof of Onsager's conjecture[END_REF], its application to fluids is still debated. In fact, the developement of finite-time singularities in Euler flows is until now an unsolved problem [START_REF] John D Gibbon | The three-dimensional euler equations: Where do we stand?[END_REF], while the same question formulated for the Navier-Stokes equations is among the open Millennium Prize Problems of the Clay Mathematics Institute [START_REF] Fefferman | Existence and smoothness of the Navier-Stokes equation[END_REF]. This debate concerns the existence of singularities in real space. In 1981, Frisch and Morf [START_REF] Frisch | Intermittency in nonlinear dynamics and singularities at complex times[END_REF] paved the way to another possibility based on the existence of complex singularities. They proved on a simple one-dimensional non-linear Langevin system that the dynamics of such complex singularities could be directly connected to intermittency, as dissipation bursts occur whenever a complex singularity approaches the real axis.

Since then, this scenario was also confirmed in the one-dimensional Burgers equation -a 1D surrogate of the Navier-Stokes equation. In this system, real singularities can be observed in the inviscid limit and manifest as shocks, i.e. finite jumps in the velocity. Shocks dissipate energy in agreement with the dissipation anomaly [START_REF] Eyink | Turbulence Theory[END_REF]. They correspond to the collapse of two complex conjugate singularities onto the real axis [START_REF] Fournier | L'équation de burgers déterministe et statistique[END_REF][START_REF] Senouf | Pole dynamics and oscillations for the complex burgers equation in the small-dispersion limit[END_REF]. When a viscosity ν is added, the singularities are repelled from the real axis, the closest one being constantly at a distance greater than O(ν 3/4 ) to the real axis. The complex singularities follow Calogero-Moser (CM) dynamics [START_REF] Calogero | Motion of poles and zeros of special solutions of nonlinear and linear partial differential equations and related solvable many-body problems[END_REF], with long-range interactions (decaying in 1/r). There is an exact mapping between such CM dynamics and the solution of the PDE, which can be described exactly via pole decomposition coupled to the integration of the CM equations [START_REF] Senouf | Pole dynamics and oscillations for the complex burgers equation in the small-dispersion limit[END_REF].

The generalization of these findings to 3D is challenging [START_REF] Frisch | Singularities of euler flow? not out of the blue![END_REF]. The computational burden to resolve the Navier-Stokes equation for a fluid with typical velocity U and length L scales like Re 3 , where Re ∼ U L/ν is the Reynolds number. Most of the earlier attempts to track complex singularities in the inviscid limit were performed using the "singularity strip" method [START_REF] Sulem | Tracing complex singularities with spectral methods[END_REF], which is based on the observation that the behaviour of the energy spectrum at large wavenumber k is dominated by the position of the singularity closest to the real axis, and decays like exp (-2δk), where δ is the imaginary part of corresponding singularity. Fitting the large wavenumber tail of the energy spectrum as a function of time, one then gets an estimate of δ(t), and a real singularity occurs when δ(t) = 0. So far, studies have only identified exponentially decaying regimes for δ(t) [START_REF] Siegel | Calculation of complex singular solutions to the 3d incompressible euler equations[END_REF] which suggests the absence of finite time blow-up. However, we cannot guarantee that this extrapolation is correct due to numerical limitations.

New perspectives on these issues were opened recently by Campolina and Mailybaev [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF], exploring fluid dynamics on log-lattices. This technique may be viewed as a generalization of the so-called "shell models" [START_REF] Gloaguen | A scalar model for mhd turbulence[END_REF][START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF] and solves the equations of motion in Fourier space using a sparse set of Fourier modes. The modes are evenly spaced points in log space ("logarithmic lattices"). They interact via nonlinear equations derived from the fluid equations by substituting for the convolution product a new operator, which can be seen as a convolution on the log-lattice, while preserving most symmetries of the original equation. The model is valid for all dimensions. In 1D, it was shown to encompass [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF] the dyaic and Sabra shell models of turbulence [START_REF] Gloaguen | A scalar model for mhd turbulence[END_REF][START_REF] Biferale | Shell models of energy cascade in turbulence[END_REF]. In 3D, its solutions have the same behaviour as the Navier-Stokes equation in Fourier space (energy spectrum, energy transfers), over an unprecedented wide range of scales [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF]. In the inviscid equations, a finite-time blow-up is observed [START_REF] Campolina | Chaotic blowup in the 3d incompressible euler equations on a logarithmic lattice[END_REF] in connection with a chaotic attractor that propagates at a constant average speed in a renormalized Fourier space, like a wave. However, Campolina and Mailybaev did not attempt to track possible complex singularities in connection with such a blow-up. This is the purpose of the present paper. In the first part, we validate the close connection between fluid dynamics on log-lattice and real fluid dynamics by focusing on the 1D Burgers equation, where dominant complex singularities are tracked using the singularity strip method. In the second part, we extend this technique to 3D to obtain new scalings regarding the approach to the real axis and the influence of normal, hypo and hyper dissipation.

Log-lattice framework

Definitions and notations

We consider a d-dimensional complex vector field u(t, k) = (u 1 , . . . , u d ) depending on time t ∈ R and on the wave vector k = (k 1 , . . . , k d ). We shall interpret u as the Fourier components of the velocity field. For this reason, we require them to satisfy the Hermitian symmetry u(t, -k) = u(t, k) with respect to k, which is the Fourier property of a real-valued function in physical space. The wave vector k is embedded on a logarithmic lattice (in short, log-lattice), which means that its components follow geometric progressions k = k 0 (±λ m 1 , . . . , ±λ m d ) for integers m 1 , . . . , m d , where k 0 = 2π is a fixed positive reference scale, and λ > 1 is the spacing factor of the lattice. The dependence of u on t and k is henceforth implicit and specified only when ambiguity prevails.

Fluid dynamics on log-lattice [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF] is the set of vector fields u which are solutions of the equations

k β u β = 0, (1a)
∂ t u α + ik β (u α * u β ) = -ik α p -νk 2γ u α + f α , (1b)
(u α * u β )(k) = q+r=k u α (q)u β (r), ( 1c
)
where p is the complex pressure field that enforces incompressibility (1a), f is a vectorial forcing, and ν is a non-negative viscosity parameter. When ν > 0, the exponent γ measures the dissipation degree: we say the flow has viscous (or usual

) dissipation if γ = 1, it has hypo-dissipation if γ < 1, and it has hyper-dissipation if γ > 1.
Similarly to the dynamics of continuous media, system (1) is the incompressible Navier-Stokes equations on the log-lattice. When ν = 0, the flow is inviscid, and the system reduces to the incompressible Euler equations on the log-lattice.

The convolution in eq. (1c) defines triadic interactions on the logarithmic lattice, which are nontrivial only if the equation λ m = ±λ q ± λ r has integer solutions m, q, r. As shown in [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF], this is possible only for particular values of λ, which determine the number of possible interactions on the grid. In this paper, we consider the following three values: λ = 2, with 3 interactions per direction; λ = φ ≈ 1.618 (the golden number), with 6 interactions per direction; and λ = σ ≈ 1.325 (the plastic number), with 12 interactions per direction. As λ decreases from 2 to σ, the density of nodes and the number of interactions on the grid increase. We recall, however, that the interactions for these log-lattices are all local.

Global quantities

By analogy with the Fourier representation of classical fluid flows, we define the global quantities representing the total energy E and the helicity H as

E = k |u| 2 , (2)
H = k u α ω α , (3)
where ω α = αβγ ik β u γ is the vorticity field; here, αβγ is the Levi-Civita symbol. Regular solutions of the unforced three-dimensional inviscid system (1) conserve these quantities in time [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF]. Moreover, we define the energy spectrum E(k) as

E(k) = |u| 2 S k , (4)
where the average • S k is taken over the wave vectors in the shell S k delimited by spheres of radii k and λk. More explicitly,

|u| 2 S k = 1 N k (λk -k) k≤|q|<λk |u(q)| 2 , (5)
where

N k ∼ (log k) d-1
is the number of wave vectors in the shell S k .

Regularity

The solutions of fluid dynamics equations on log-lattices (1) share some regularity properties with the original models. The main mathematical results are for the inviscid Euler equations [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF]. For this system, the local-in-time existence of strong solutions and a Beale-Kato-Majda blow-up criterion were proved. Exploiting the conservation of enstrophy, one proves the global regularity of two-dimensional flows. In the threedimensional case, high-resolution numerical log-lattice simulations disclosed a finite-time blow-up, characterized by a chaotic wave travelling with constant average speed along a renormalized set of variables [START_REF] Campolina | Chaotic blowup in the 3d incompressible euler equations on a logarithmic lattice[END_REF]. Such blow-up scenario was confirmed for λ = φ and λ = σ, presenting the same asymptotic blow-up scalings [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF]. In the viscous case, numerical simulations suggest the expected global regularity of solutions.

Singularity strip method for log-lattices

If a potential singularity is due to an imaginary pole crossing the real axis, one can track its distance to the real axis via the singularity strip method [START_REF] Sulem | Tracing complex singularities with spectral methods[END_REF]. This method considers the analytic continuation u(z) of the physical-space velocity field and is based on the following property: if

u(z) ∼ 1/(z -z * ) ξ , for z → z * (6)
in a neighborhood of the complex singularity z * = a + iδ, then its Fourier transform ûk satisfies ûk ∼ k -d-ξ e ika e -δk , as k → ∞.

Asymptotics of [START_REF] Fournier | L'équation de burgers déterministe et statistique[END_REF] provide the corresponding exponential decay E(k) ∼ e -2δk for the energy spectrum over a typical length 2δ. Therefore, one can measure the distance of the dominant pole to the real axis by monitoring the decay of the energy spectrum in Fourier space. A finite-time singularity at instant t b would occur if δ → 0 as t → t b . Extension of this notion to the log-lattice framework is natural. It relies on the observation that if a flow (1) on log-lattice satisfies u(k) ∼ k -d-ξ e -δk , then its inverse Fourier transform obeys a relation similar to [START_REF] Eyink | Turbulence Theory[END_REF]. Therefore, we can generalize the singularity strip method to log-lattices, where 2δ is estimated from the slope of log E(k) as a function of k.

Numerical methods

Equations (1) are numerically integrated using a technique analogous to viscous splitting. Considering a time step dt, we obtain u(t + dt) from u(t) employing the following strategy. Using u(t) as initial condition, we first solve the inviscid equation

∂ t u α = P αβ [-ik σ (u β * u σ ) + f β ] , (8)
where

P αβ = δ αβ - kαk β k 2
accounts for the pressure term under the incompressibility hypothesis (1a). For that, we use an explicit 4th order Runge-Kutta method. This yields u(t + dt) ν=0 . Then, the viscosity is taken into account through u(t + dt) = u(t + dt) ν=0 e -νk 2γ dt . In this whole process, we adapt dynamically the time step dt.

1D Burgers equation

Before going to the full three-dimensional Navier-Stokes system on log-lattices, we take an intermediate step by studying the simpler one-dimensional Burgers equation. For this system, several exact mathematical results are available. This allows us to probe the singularity strip method on log-lattices, by comparing our numerical computations with the exact expected results.

The one-dimensional Burgers equation on log-lattices is obtained from system (1) as follows. We consider a compressible pressureless flow on a one-dimensional log-lattice. Mathematically, this translates into setting p = 0 and dropping eq. (1a) from the system, which reduces to

∂ t u + u * ∂ x u = -νk 2γ u + f, (9a)
(u * ∂ x u)(k) = q+r=k ir u(q)u(r). (9b)
It was shown [START_REF] Campolina | Fluid Dynamics on Logarithmic Lattices and Singularities of Euler Flow[END_REF] that, up to a prefactor in the convolution (9b), the Burgers equation on log-lattices is equivalent to well-known shell models of turbulence for specific choices of parameters. Particularly, when λ = 2, system (9) (but with a factor 2 added in the convolution and restricting to imaginary solutions) is the dyadic model [START_REF] Desnyansky | The evolution of turbulence spectra to the similarity regime[END_REF], while for λ = φ (but with a factor -φ 2 added in the convolution) it is the Sabra model [START_REF] L'vov | Improved shell model of turbulence[END_REF] in a three-dimensional parameter regime (second invariant is not sign defined). Because of this relation with shell models of turbulence, the Burgers equation on the onedimensional log-lattice inherits several results concerning the regularity of its solutions, which we briefly review now.

For the dyadic model (λ = 2) with ν > 0, there are theorems [START_REF] Cheskidov | Blow-up in finite time for the dyadic model of the navier-stokes equations[END_REF] for global existence of weak solutions (satisfying the energy inequality at almost all time), local regularity when γ > 1/3, global regularity when γ ≥ 1/2, and finite-time blow-up when γ < 1/3 for sufficiently large initial conditions. Note that for the continuous version of the 1D Burgers equations, global existence and analycity holds whenever γ ≥ 1/2, while finite-time blow-ups are present whenever γ < 1/2 [START_REF] Kiselev | Blow up and regularity for fractal Burgers equation[END_REF]. In contrast, we have presently no rigorous statements about the dyadic model for the parameter range 1/3 ≤ γ < 1/2. This means that the mathematical techniques used in the currently available theorems are not sharp enough to separate the finite-time blow-up and the global regularity regimes. The finite-time singularity in the inviscid case was also rigorously established [START_REF] Katz | Finite time blow-up for a dyadic model of the euler equations[END_REF].

For the viscous Sabra model (λ = φ) with usual dissipation γ = 1, there are proofs [START_REF] Constantin | Analytic study of shell models of turbulence[END_REF] of global regularity of strong solutions. Like the Navier-Stokes equations, the dynamics of the Sabra model develops within finite degrees of freedom. Indeed, the finite dimensionality of the global attractor and the existence of a finite-dimensional Table 1: Exponents of the inviscid scalings of various quantities measured for the 1D Burgers and the 3D Euler equations within different values of the grid spacing λ. The scalings are with respect to τ = 1-t/t b , where t b is the blow-up time. By definition, the energy spectrum scales like E(k) ∼ k -1-2α , the maximum value of the vorticity scales like ω max ∼ τ -β , and the width of the singularity strip scales like δ ∼ τ µ . The (0) inertial manifold were proved [START_REF] Constantin | Analytic study of shell models of turbulence[END_REF]. On the other hand, the inviscid model has [START_REF] Constantin | Regularity of inviscid shell models of turbulence[END_REF] global-in-time existence of weak solutions with finite energy, local-in-time regularity, and a Beale-Kato-Majda blow-up criterion. Despite the absence of rigorous proofs, it is well-known [START_REF] Mailybaev | Spontaneously stochastic solutions in one-dimensional inviscid systems[END_REF] that Sabra (in the three-dimensional parameter regime) develops a selfsimilar finite-time blow-up, characterized as a travelling wave in a renormalized system of variables (cf. [START_REF] Dombre | Intermittency, chaos and singular fluctuations in the mixed Obukhov-Novikov shell model of turbulence[END_REF]). Following the dynamical systems approach, such blow-up can be seen as a fixed-point attractor of the associated Poincaré map [START_REF] Mailybaev | Bifurcations of blowup in inviscid shell models of convective turbulence[END_REF].

To our knowledge, there are no systematic results about the development of singularities in Sabra with general dissipation exponents γ, nor in the case of our third lattice parameter λ = σ.

Inviscid flow

We start with the inviscid (ν = 0) Burgers equation with and without forcing. When forcing, initial conditions are equal to zero, and the forcing is equal to the imaginary unit i on the first mode (k 0 ) for λ = 2, the first two modes when λ = φ, and the first three modes when λ = σ. Without forcing, initial conditions are taken such that total energy E = 1, and first mode, two first modes or three first modes have positive uniform real value, depending on λ = 2, φ or σ, while initial smaller scales are zero. We observe finite-time blow-up for all three values of λ in the two cases. Numerical results are plotted in fig. 1, and scaling exponents are summarized in table 1.

The maximum of the gradient ω max (t) = max k |ku(k)| blows up in finite time, following the self-similar law

t b ω max ∼ 1 τ , τ = 1 - t t b , (10)
displayed in fig. 1b. While the blow-up time depends on the forcing and the value of λ, the self-similar law (10) is independent of these variables. This law also holds for the original continuous model. Indeed, differentiating the classical Burgers equation ∂ t u + u∂ x u = 0 with respect to x, we get that the space derivative ω = -∂ x u obeys dω/dt = ∂ t ω + u∂ x ω = ω 2 , whose solution is exactly eq. ( 10) with t b = 1/ω(t = 0). We also check that the energy spectrum evolution is universal, in the sense that it only depends on τ . This is illustrated in fig. 1a, where spectra for different λ but similar τ are shown to overlap. As τ approaches zero, the energy spectrum gradually widens towards larger values of k, developing a power-law E(k) ∼ k -2α-1 with 1 + 2α = 1.733, which corresponds to the scaling law u(k) ∼ k -α with α = 0.367. Such asymptotics agrees with exact results from the renormalization group formalism applied to the Sabra shell model [START_REF] Fontaine | Functional renormalisation group approach to shell models of turbulence[END_REF].

Finally, we compute the analyticity strip width δ as the solutions approach the blow-up. This is done using the formula [START_REF] Fournier | L'équation de burgers déterministe et statistique[END_REF] with ξ + 1 = α. The result is shown in fig. 1c. We verify that δ decays to zero in finite time, following a power law δ ∼ τ µ , with µ = 1.546. This decay is also universal and does not depend on the value of λ or the forcing. The width of the analyticity strip is closely associated with k max , defined as the wavenumber at which ω attains its maximum value. Indeed, we see in fig. 1d that k max δ is approximately constant in time. This is in agreement with the asymptotic eq. ( 7), which implies that ω max is achieved at k max ∼ 1/δ.

The self-similar law [START_REF] Frisch | Singularities of euler flow? not out of the blue![END_REF] is valid for all λ in average only. The figures show that the blow-up looks truly self-similar only for the values λ = 2 and λ = φ. The oscillations in the case λ = σ suggest a different blow-up scenario (e.g. quasi-periodic or chaotic). A detailed analysis of this is left for future work.

Viscous dissipation

We now introduce viscous dissipation (γ = 1) and study how the dynamical behaviour depends upon the viscosity parameter ν. This section restricts the analysis to the value λ = 2. We introduce a force at the large scale, whose amplitude is adapted dynamically so that the total power input is constant in time (f k=k 0 = P u k=k 0 /|u| 2 k=k 0 , where P = 1). In this setup, the dissipative term is strong enough to prevent the blow-up, and the solution reaches stationarity. The energy spectrum develops a power law in the intermediate scales (called the inertial range) followed by an exponential decay at larger k -see fig. 2a. In the inertial range,

E(k) ∝ k -5/3 , corresponding to u(k) ∝ k -1/3 .
The maximum value of the derivative ω max is inversely proportional to the viscosity, following the power law ω max ∼ ν -β with β = 0.5, as shown in fig. 2b. This scaling law can be derived when assuming a viscosity-independent anomalous dissipation > 0 in the inviscid limit ν → 0. Under this assumption, we have the balance νω 2 ∼ , which provides ω ∼ ( /ν) 1/2 .

Accordingly, the width of the analyticity strip does not decline to zero. However, it stabilizes at a finite value that depends on the viscosity -see fig. 2c -and follows the power-law scaling δ ∼ ν µ , with exponent µ = 0.7067. This is smaller than expected from a dimensional argument "a la Kolmogorov", in which = νu 2 /δ 2 , with u ∼ δ 1/3 , would instead predict δ ∼ ν 3/4 . The strip width follows approximately the scaling δ ∼ 1/k max , as shown in fig. 2d.

Hyper-and hypo-dissipation

We have also studied the influence of the dissipation degree γ on the various scaling laws. This is summarized in fig. 2 andtable 2. The slope of the spectrum is insensitive to γ and displays a E(k) ∼ k -5/3 law with no intermittency correction. On the other hand, the slopes of both the inverse of the maximum gradient and the singularity width increase in absolute value as γ is decreased towards 1/3. We defer the discussion about those results to section 5.1.

Table 2: Scaling exponents of various quantities as a function of γ measured for the 1D Burgers and the 3D Navier-Stokes equations with grid spacing λ = 2. The scalings are with respect to the viscosity ν. By definition, the energy spectrum scales like E(k) ∼ -1-2α , the maximum value of the vorticity scales like ω max ∼ ν -β and the width of the analyticity strip scales like δ ∼ ν µ .

1D Burgers

3D Navier-Stokes

γ α β µ α β µ 1/3 - - - 2

Critical dissipation degree γ = 1/3

According to [START_REF] Cheskidov | Blow-up in finite time for the dyadic model of the navier-stokes equations[END_REF], there are finite time blow-up solutions for the Burgers equation ( 9) with λ = 2 whenever γ < 1/3. However, the theorems say nothing about the limit case γ = 1/3. For this reason, we call this value as being the critical dissipation degree. It is natural to ask whether the blow-up might or might not occur in this specific situation.

Here we consider not only λ = 2, but also extend this question to the other two lattice parameters.

We initialized the flow with the same data as in the inviscid case and set the small viscosity ν = 10 -7 . We observed a finite time blow-up for all three λ, illustrated in fig. 3. The blow-up time is larger than in the inviscid case, but the scaling laws are the same -both the prefactor and the scaling exponents -as in the inviscid case. The exponents are summarized in table 3. The only exception is for the scaling law of δ in the case λ = σ. This might be due to the oscillations in the energy spectrum, making it harder to fit the exponential decreasing, see fig. 3a. In this small-viscosity run, viscosity only delays the blow-up but does not influence the development of the singularity. However, we observed a surprising behaviour change when increasing the viscosity to larger values. There is a transition between a smallviscosity regime, where finite time blow-up occurs, and a large-viscosity regime, where the blow-up disappears, and the width of the analyticity strip saturates to a finite value -see insert of fig. 3c. The amplitude of δ seems to follow a critical mean-field behaviour, as it varies like δ ∼ √ ν -ν c , with ν c ∼ 0.4. A similar transition is observed at a lower value of γ, with ν c increasing as γ decreases.

This transition is in fact not contradicting the mathematical results by [START_REF] Cheskidov | Blow-up in finite time for the dyadic model of the navier-stokes equations[END_REF], since they prove existence of blow-up only for initial conditions larger than a threshold that depends linearly in the viscosity. In all our calculations, we start with the same initial conditions. This means that for large enough values of viscosity, the initial condition becomes smaller than the threshold, therefore invalidating the hypothesis of the theorem. More than that, our numerical results suggest that this hypothesis is actually essential for the result of the theorem and might not be dropped in general.

3D Euler and Navier-Stokes equations

Inviscid flow -Euler equations

We now turn to the full three-dimensional incompressible fluid dynamics on log-lattices, starting with the inviscid Euler equations. We consider here the three lattice spacings λ. In order to test universality, we ran the case λ = φ with two different incompressible random initial conditions, differing by their range of scales. Default initial conditions are defined at large scale |k| < 3k 0 , while the other (denoted by a star * ) are defined at scales |k| < k 0 λ 3 . We observed a finite-time blow-up in all setups, in agreement with previous results documented in [START_REF] Campolina | Chaotic blowup in the 3d incompressible euler equations on a logarithmic lattice[END_REF][START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF]. Here, we observe that while the blow-up time depends on the initial conditions, the dynamics become universal when plotted in non-dimensional variables, as illustrated in fig. 4. The spectra for distinct values of λ overlap when plotted at the same non-dimensional times τ = 1 -t/t b , as evidenced in fig. 4a. The slope of the power law in the inertial range is steeper than in 1D Burgers, with a value very close to -7/3. This is the slope expected for a helicity cascade. Our exponent is slightly smaller than those found in some direct numerical simulations of the Euler equations, where a E(k) ∼ k -3 spectrum is observed [START_REF] Orlandi | Vortex events in Euler and Navier-Stokes simulations with smooth initial conditions[END_REF][START_REF] Bustamante | Interplay between the beale-kato-majda theorem and the analyticity-strip method to investigate numerically the incompressible euler singularity problem[END_REF], but comparable to the value 2.33 obtained in more recent simulations [START_REF] Fehn | Numerical evidence of anomalous energy dissipation in incompressible Euler flows: towards grid-converged results for the inviscid Taylor-Green problem[END_REF].

The maximum value of the vorticity ω max diverges during the blow-up, as shown in fig. 4b. Its asymptotic scaling is the same as for the maximum gradient in the 1D Burgers equation, given by eq. [START_REF] Frisch | Singularities of euler flow? not out of the blue![END_REF]. However, contrarily to the 1D case, the constant in front of the power law varies as a function of λ and is not simply given by 1/t b . This is not too surprising given the 3D nature of the flow, which prevents the application of the simple blow-up argument used for 1D Burgers. However, as λ is decreased towards 1, the non-dimensional curve becomes closer to the exact asymptotic law. Approaching the blow-up, the width of the analyticity strip decays to zero with a power law δ ∼ τ µ with exponent µ ≈ 2.81 -see fig. 4c. This is larger than in 1D Burgers. This decay is also universal, as it does not depend on λ. However, it does not show a simple dependence with k max as seen in fig. 4d. This might be related to the chaotic nature of the blow-up attractor [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF].

Viscous dissipation -Navier-Stokes equations

We now add the viscous term with γ = 1 and a constant-power forcing. The solutions achieve a statistically stationary state in this framework, whose average scalings are depicted in fig. 5. Like in 1D Burgers, the energy spectra display a power law until the solution reaches the viscous scale, with an inertial range widening as ν decreases. The slope of the energy spectrum is slightly steeper than Kolmogorov's -5/3, with an intermittency correction of around 0.13. Accordingly, the scaling exponent α for u(k) ∼ k -α is α = 0.40, slightly larger than 1/3. The maximum vorticity ω max increases with decreasing viscosity, following the power law ω max ∼ ν -β with an exponent β = 0.39 lower than Kolmogorov's 1/2.

The width of the analyticity strip decays with viscosity as δ ∼ ν µ with an exponent µ = 0.65 -see fig. 5c. Such decay is less intense than in the 1D Burgers equation. Nevertheless, the dependence of δ on 1/k max in the Navier-Stokes case is sharper, as one verifies by comparing fig. 5d with fig. 2d.

Hyperviscous dissipation

We now consider what happens in the hyperviscous case γ > 1. We keep the constantpower forcing to reach stationary states.

For γ = 2, we still observe a power-law energy spectrum followed by an exponential cut-off at the viscous scales -see fig. 5a. The inertial range keeps widening as ν is decreased. The slope of the energy spectrum is very close, but slightly steeper than -5/3. The exact fitting provides us an intermittency correction around 0.07, corresponding to α = 0.37, see table 2. The maximum vorticity ω max increases with decreasing viscosity like a power law, with an exponent β = 0.19 lower than usual (γ = 1) viscous case. The width of the analyticity strip decays with viscosity with an exponent µ = 0.26. Like in the viscous case, δ appears to scale simply like 1/k max , as seen on fig. 5d.

The above results suggest that the intermittency corrections in the energy spectra are smaller for hyperdissipation. Indeed, as the dissipation degree γ increases, the exponent α converges towards Kolmogorov's 1/3, see table 2. We checked that for the stronger degree γ = 8, they vanish completely, and the dependence of ω max and δ on ν become very weak. This is explained by the very sharp viscous cut-off due to the hyperviscous dissipation. Indeed, the equivalent of the Kolmogorov scale k d in the hyperviscous case relates to ν as k d ∼ ν 1/(1-1/3-2γ) , becoming independent of viscosity in the limit γ → ∞. For γ = 2 the dependence is δ ∼ k -1 d ∼ ν 0.3 , close to what is observed for the scaling of the singularity strip width.

Hypoviscous dissipation

The case with hypoviscous dissipation 1/3 < γ < 1 is qualitatively similar to the viscous and hypervisous cases -see fig. 5. Exponents, however, are steeper. The corresponding values are reported in table 2. The energy spectrum develops a slope corresponding to the exponent α = 0.5, which is steeper than Kolmogorov's 1/3 but milder than Euler's 2/3 on log-lattices. The singularity width appears again to be controlled by the wave number corresponding to the maximum vorticity -see fig. 5d. On the other hand, the maximum vorticity grows much more rapidly than in the viscous case, with an exponent twice as big, as shown in fig. 5b. This may indicate that we are approaching a critical dissipation degree, below which finite-time blow-up will occur.

Critical dissipation degree γ = 1/3

The asymptotics of Kolmogorov's length scale for a flow with a general dissipation degree predicts the breakdown of the viscous cut-off when γ approaches the critical value 1/3. Indeed, the dissipation scale k d is obtained from the dimensional balance between the convective and the dissipative terms k d u 2 d ∼ νk 2γ d u d . On the other hand, Kolmogorov's theory states that u d ∼ 1/3 k -1/3 d for the energy dissipation rate , which has a finite positive value in the inviscid limit. Together, these expressions yield

k d ∼ 1 6(γ-1/3) ν 1 2(1/3-γ) , (11)
which, for sufficiently small ν, provides k d → +∞ when γ 1/3. For this reason, we call γ = 1/3 the critical dissipation exponent, the value at which we expect that the dissipative term is no longer strong enough to prevent a finite-time singularity. We recall this was the case for the 1D Burgers equation on log-lattices.

Motivated by the above arguments, we investigate the critical hypo-diffusive degree in the full 3D system on log-lattices. The following analysis considers the spacings λ = 2 and λ = φ. The initial data is the same as we used in the inviscid simulations, and viscosity is the same ν = 10 -7 .

In this regime, we observed a finite time blow-up for the two values of λ, illustrated in fig. 6. Like in 1D Burgers, the blow-up time is larger than in the inviscid case, but the scaling laws are the same. This is summarized in table 3. The slope of the energy spectrum remains -7/3. For λ = 2 and ν = 10 -3 , the dynamics becomes stationary, meaning there is as in 1D a phase transition, but between ν = 10 -3 and 10 -7 , smaller than ν c ∼ 0.4 in 1D.

Discussion

Scaling laws

The variations of the scaling exponents with respect to the diffusion exponent γ are shown in fig. 7.

Predictions for the scaling laws are possible using simple dimensional arguments if we impose δ ∼ 1/k max , as empirically observed. Indeed, from u ∼ k -α and ω ∼ ku, we get ω max ∼ k 1-α max ∼ δ α-1 so that we get:

β = µ(1 -α). (12)
This fixes a link between the 3 exponents that is well satisfied -see fig. 7b. On the other hand, one can connect µ and α by extending the argument fixing the Kolmogorov scale to hypo and hyper-viscous cases: we impose that k max is fixed by the condition that the viscous term balances the non-linear term νk 2γ max u max ∼ k max u 2 max . Using u max ∼ k -α max and δ ∼ 1/k max we then get: This prediction is tested in fig. 7c and is well satisfied. Without loss of generality, the only free parameter can be taken as α(γ). In the limit γ → 1/3, we can fix it by imposing that β = 1, which is the scaling corresponding to conservation of the circulation of u [START_REF] Pomeau | On the self-similar solution to the euler equations for an incompressible fluid in three dimensions[END_REF]. From eqs. ( 12) and ( 13), we then get α = 1 -γ = 2/3, corresponding to a helicity cascade. In all other cases, we have no clear theories to predict the variations of α with γ. Notably, when γ → ∞, we recover α = 1/3 corresponding to an energy cascade. . The data points are reported from tables 2 and 3, while the dotted lines correspond to eq. ( 12); (7c) For the scaling of the singularity strip δ ∼ ν µ . The data points are reported from tables 2 and 3. while the dotted lines correspond to eq. ( 13).

µ = - 1 1 -α -2γ , (13)

Interest of the critical case

The critical case γ = 1/3 is more than purely academic: renormalization group (RNG) analysis of NSE in Fourier space [START_REF] Yakhot | Renormalization-group analysis of turbulence[END_REF] indeed shows that the fixed point of the equations corresponds to a Navier-Stokes equation with turbulent viscosity scaling like A 1/3 k -4/3 , where A is a constant with value A = 0.1447 in 1D and A = 0.4926 in 3D. This corresponds exactly to eq. ( 1), with γ = 1/3 and ν = A 1/3 . This model is sometimes used as a subgrid model of turbulence [START_REF] Laval | Nonlocality and intermittency in three-dimensional turbulence[END_REF]. In that respect, it is interesting that the transition viscosity found in sections 3.4 and 4.5 (at constant injected power, i.e. = 1) is very close to the RNG value in 1D. On the one hand, this guarantees that the size of the inertial range is very wide, in agreement with the RNG picture of scale invariant solutions. On the other hand, this means that the solution is very close to a blow-up, which could have implications regarding the stability of this subgrid scheme.

Implications for real Euler or Navier-Stokes?

The log-lattices simulations we performed cannot be seen as an exact model of the Euler or Navier-Stokes equations because they remove by construction many non-linear interactions of the original equations, especially the non-local one. Nevertheless, because they obey the same conservation laws and symmetries, they may capture some scaling laws of the original equation more accurately. Comparing our findings with the few results on the topic is engaging.

Regarding the Euler equation, recent high-resolution numerical simulation in the axisymmetric case by [START_REF] Swetha | Insights from a pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible euler equation[END_REF] explored the scaling of the singularity strip in the blowing situation proposed by [START_REF] Luo | Potentially singular solutions of the 3d axisymmetric euler equations[END_REF]. They found an exponent µ = 2.6 ± 0.5, which is compatible with the value 2.8 ± 0.1 that we get from table 2. Unfortunately, they do not provide an estimate of the slope of the energy spectrum. Previous older results in the Taylor-Green vortex [START_REF] Brachet | Numerical evidence of smooth self?similar dynamics and possibility of subsequent collapse for three?dimensional ideal flows[END_REF][START_REF] Bustamante | Interplay between the beale-kato-majda theorem and the analyticity-strip method to investigate numerically the incompressible euler singularity problem[END_REF] found a steeper spectrum corresponding α ∼ 1. However, spectra with exponent matching our -7/3 value were observed in the early stage of recent simulations at larger resolution [START_REF] Fehn | Numerical evidence of anomalous energy dissipation in incompressible Euler flows: towards grid-converged results for the inviscid Taylor-Green problem[END_REF]. Therefore, the main characteristics of blow-up in log-lattices simulations agree with the most recent results observed in the traditional DNS of the Euler equation.

Regarding the Navier-Stokes equations, we can look at two recent results. The first one by [START_REF] Dubrulle | Multi-fractality, universality and singularity in turbulence[END_REF] finds a value of µ = 0.89 using recent DNS of NSE. This value is larger than the value we found in the present paper, corresponding to µ = 0.65. Another recent result [START_REF] Buaria | Extreme velocity gradients in turbulent flows[END_REF] estimates β in 3D NSE. They indeed found that the tail of the PDF of enstrophy scales like ν 0.77 τ -2 K , where τ K ∼ ν 1/2 is the Kolmogorov time. Identifying such extreme events of enstrophy with ω 2 max , we thus get β DN S ∼ 0.88, which is also much larger than the value we observe in log-lattices β LL ∼ 0.39. Note, however, that both DNS values are compatible with eqs. ( 12) and ( 13), provided we choose α ∼ 0, hinting at the presence of multifractality. Log-lattices simulations are generally much less intermittent than DNS [START_REF] Ciro | Fluid dynamics on logarithmic lattices[END_REF], with one dominating exponent (monofractal behaviour). Some time ago, [START_REF] Laval | Nonlocality and intermittency in three-dimensional turbulence[END_REF] linked the intermittency properties of NSE with non-local interactions, which is coherent with this observation. Therefore, the difference between log-lattices simulations and DNS could be explained by differences in the amount of non-local interactions.

Figure 1 :

 1 Figure 1: Inviscid blow-up for the 1D Burgers equation for λ = 2 (yellow), λ = φ (red) and λ = σ (blue). Continuous lines and filled symbols indicate simulations with constant forcing, while dotted lines and open symbols indicate simulations without forcing. (1a) Spectra at different renormalized relative time τ = 1 -t/t b . The black dotted line has a slope of -1.733. (1b) Maximum value of the derivative 1/t b ω max as a function of τ ; The black dotted line is the theoretical value from eq. (10). (1c) Width of the analyticity strip 2δ as a function of τ . The black dotted line has a slope given in table 1. (1d) Renormalized width k max δ as a function of τ . The black dotted line has a slope of 0.

Figure 2 :

 2 Figure 2: Influence of the type of viscosity on the stationary dynamics of the viscous 1D Burgers equation for λ = 2 and γ = 0.5 (hypo-viscous case, blue circle), γ = 1 (viscous case, red squares) and γ = 2 (hyperviscous case, yellow diamond). (2a) Energy spectrum. The black dotted line has a slope -5/3; (2b) Maximum value of the derivative 1/ω max as a function of viscosity. The black dotted line has a slope given in table 2 for each case. (2c) Width of the analyticity strip 2δ as a function of viscosity. The black dotted line has a slope given in table 2 for each case. (2d) Renormalized width k max δ as a function of viscosity.

Figure 3 :

 3 Figure 3: Blow-up for the critical (γ = 1/3) 1D Burgers equation, with ν = 10 -7 and for λ = 2 (yellow), λ = φ (red) and λ = σ (blue). (3a) Spectra at different renormalized relative time τ = 1-t/t b ;(3b) Maximum value of the derivative 1/t b ω max as a function of τ (3c) Width of the analyticity strip 2δ as a function of τ . The insert shows the behaviour of the width of the analyticity strip at t = ∞ when the viscosity is increased, for γ = 1/3 (blue data points) (resp. γ = 1/4 (red data points)). The dotted lines are fits of the type √ ν -ν c , with ν c = 0.4 (resp. ν = 0.9). (3d) Renormalized width k max δ as a function of τ . The dotted line has the same scaling and prefactor as in the inviscid blow-up case, see fig.1

Figure 4 :

 4 Figure 4: Inviscid blow-up for the 3D Euler equations for λ = 2 (yellow), λ = φ (red) and λ = σ (blue). (4a) Spectra at different renormalized relative time τ = 1 -t/t b , from 0.2542 to 0.00001 from left to right. Spectra with continuous lines and dotted lines correspond to different initial conditions. The black dotted line has a slope of -7/3; (4b) Maximum value of derivative 1/t b ω max as a function of τ . The black dotted line has a slope of 1; (4c) Width of analyticity strip 2δ as a function of τ ; The black dotted line has a slope 2.805. (4d) Renormalized width k max δ as a function of τ . In panels 4b, 4c and 4d, we used different symbols for different initial conditions: circles, and squares.

5 :

 5 Stationary dynamics for the 3D Navier-Stokes equations for λ = 2 and γ = 0.5 (hypo-viscous case, blue circle), γ = 1 (viscous case, red squares) and γ = 2 (hyperviscous case, yellow diamond). (5a) Energy spectrum. The black dotted line has a slope -5/3; (5b) Maximum value of the derivative 1/ω max as a function of viscosity. The black dotted line has a slope given in table 2 for each case. (5c) Width of the analyticity strip 2δ as a function of viscosity. The black dotted line has a slope given in table 2 for each case. (5d) Renormalized width k max δ as a function of viscosity.

Figure 6 :

 6 Figure 6: Blow-up for the critical (γ = 1/3) 3D Navier-Stokes equations, for λ = 2 (yellow) and λ = φ (red). (6a) Spectra at different renormalized relative time τ = 1 -t/t b ; (6b) Maximum value of the derivative 1/t b ω max as a function of τ (6c) Width of the analyticity strip 2δ as a function of τ (6d) Renormalized width k max δ as a function of τ . The dotted line has the same scaling and prefactor as in the inviscid blow-up case, see fig. 4

Figure 7 :

 7 Figure7: Variation of the scaling exponents as a function of the dissipation degree γ for the 1D Burgers (blue circle) and the 3D Navier-Stokes (red square) equations with λ = 2. (7a) For the scaling of the velocity u ∼ k -α . (7b) For the scaling of the maximum vorticity ω max ∼ ν -β . The data points are reported from tables 2 and 3, while the dotted lines correspond to eq. (12); (7c) For the scaling of the singularity strip δ ∼ ν µ . The data points are reported from tables 2 and 3. while the dotted lines correspond to eq. (13).

  superscript indicates a simulation performed with no forcing. The * superscript indicates a simulation made with a different initial condition.

		λ	t b	α	β	µ
		2	0.3898 0.37 1 1.55
		φ	0.5193 0.37 1 1.55
	1D Burgers	σ	0.4300 0.37 1 1.55
		2 (0) 0.2687 0.37 1 1.55
		φ (0) 0.1460 0.37 1 1.55
		2	0.8481 0.67 1 2.81
	3D Euler	φ φ * 0.1542 0.67 1 2.82 5.8005 0.67 1 2.83
		σ * 0.8430 0.67 1 2.67

Table 3 :

 3 Exponents

	/3	1	2.81

in the critical case γ = 1/3 of various quantities measured for the 1D Burgers and the 3D Navier-Stokes equations with different values of the grid spacing λ. The scalings are with respect to τ = 1 -t/t b , where t b is the blow-up time. By definition, the energy spectrum scales like E(k) ∼ k -1-2α , the maximum value of the vorticity scales like ω max ∼ τ -β and the width of the analyticity strip scales like δ ∼ τ µ .

Acknowledgments

This work received funding from the Ecole Polytechnique, from ANR TILT grant agreement no. ANR-20-CE30-0035, and from ANR BANG grant agreement no. ANR-22-CE30-0025. CC is thankful for the financial support from CEA during his visit to Paris-Saclay, where this work was elaborated. AM was supported by CNPq grant 308721/2021-7 and FAPERJ grants E-26/201.054/2022, E-26/210.874/2014.


Ewa Napieralska-Juszczak
email: ewa.napieralskajuszczak@univ-artois.fr

Daniel Roger

K Komeza

N Takorabet

DESIGN OF A PERMANENT MAGNET MACHINE ABLE TO WORK AT HIGH TEMPERATURES

The paper presents a methodology for assessing the effects of using different materials for magnetic wedges and applying magnet segmentation on the power losses and performance of permanent magnet machines operating at high temperatures. Time stepping finite element simulation in 2D and 3D is employed to obtain field and current distributions, power losses and derive torque ripples under various configurations.

I. INTRODUCTION

Aircrafts designed to have lower environmental footprints will use more electric power that will replace hydraulic power and hot air sources. The aircraft industry will need compact and reliable generators and actuators able to work in a very hot environment, especially near the propulsion turbines. With the classical polymer-based winding technology the internal temperatures cannot exceed about 280°C for steady state operations. In order to work continuously at higher internal temperatures it is necessary to use inorganic Electrical Insulation Systems (EIS) and to adapt the machine design to this inorganic technology [START_REF] Cozonac | Materials characterization and geometry of a high temperature induction machine[END_REF][START_REF] Hooker | High-temperature electrical insulation for egs downhole equipment[END_REF]. The inorganic insulation layers have poor mechanical characteristics and the coils must be rigid objects placed on the stator teeth. The structure of the motor winding must be a fractional-slot concentrated winding. Two structures are available for rotor: for the first one, the magnets are placed on the external surface or the rotor; for the second one, the magnets are buried in soft magnetic parts of the rotor. To ensure an adequate magnetic behaviour at high temperature (HT°), special permanent magnets were chosen characterized by relatively large electrical conductivity in comparison to standard magnets. For this reasons a special attention must be paid to the eddy currents in permanent magnet. To analyse this phenomenon the time-stepping 2D and 3D finite element method (FEM) was used. The maximum continuous operating temperature of permanent magnets depends on the temperature coefficients TC(HcJ) of the coercivity, TC(Br), of the remanent polarization, and the load line B/μ0H of the magnet. The high temperature necessitates the use of permanent magnets with a high Curie point. In our case the chosen material is Vacomax Sm 2 Co 17 [START_REF] Rodewald | Topology of Nd-F-B magnets with a high energy density[END_REF], characterised by high electrical conductivity compared with ferrites. The numbers of poles and slots were carefully chosen to reduce the high frequency slotting effects on PM losses and torque ripples. A popular method to limit the ripple torque and eddy current losses in PMs involves the use of magnetic wedges. The numbers of poles and slots were carefully chosen to reduce the high frequency slotting effects on PM losses and torque ripples [START_REF] Tang | Analysis of Iron Loss in Interior PM Machines With Distributed Windings Under Deep Field Weakening[END_REF][START_REF] Pfister | Slotted Permanent-Magnet Machines: General Analytical Model of Magnetic Fields, Torque, Eddy Currents, and Permanent-Magnet Power Losses Including the Diffusion Effect[END_REF][START_REF] Wrobel | Winding Design for Minimum Power Loss and Low-Cost Manufacture in Application to Fixed-Speed PM Generator[END_REF][START_REF] Ruoho | Temperature Dependence of Resistivity of Sintered Rare-Earth Permanent-Magnet Materials[END_REF]. A popular method to limit the ripple torque and eddy current losses in PMs involves the use of magnetic wedges. Three different designs have been considered in this paper: without magnetic wedges (NM), with magnetic wedges build from the same material as the stator core (M), and semi-magnetic wedges made of a composite material (SM).

II. SIMULATION METHOD

For the purpose of the 3D simulation a quarter of the machine has been considered to benefit from the magnetic symmetry. The windings are modelled by connections to external resistances. The time stepping analysis was performed assuming constant speed of rotation of the rotor. The 3D mesh contains 4.264778 million elements in total. In the magnets, as the conductive areas, the elements size was chosen to be less than one-third of the penetration depth. Moreover, the distribution of the density of the mesh is greater than to correctly reproduce the effect of rotation of the rotor relative to the stator. To reduce the computation times the simulations were performed on a virtual machine with eight cores First the simulation was performed for the surface magnets machine. The current density distribution on and close to the surface of the permanent magnets is shown in Fig. 1. As can be seen, the skin depth is small; hence the total losses in the magnet are not high, despite the large value of the surface current. However, the permanent magnet losses cause significant concern asdue to the design of the machinethe removal of heat from the rotor is inhibited. Excessive losses in the rotor may therefore result in local temperatures exceeding the Curie point and hence causing an irreversible demagnetization of parts of the magnet. Some limitation of the losses in the magnets may be accomplished by magnet segmentation [START_REF] Huang | Optimization of Magnet Segmentation for Reduction of Eddy-Current Losses in Permanent Magnet Synchronous Machine[END_REF]. Magnet segmentation results in significant reduction of PM losses. Design enables a significant reduction in losses in the magnets is the design of the magnets embedded into the rotor known as interior PM. Therefore, with similar stator winding, the voltage induced in the phase is smaller than in the case of surface mounted magnets. Figure 2 shows the stator voltage, the cogging torque at no-load and the comparison with the case of surfacemounted magnets. The benefits of the embedded structure are a smoother voltage and a much smaller cogging torque even than for the construction of non-magnetic wedges. For the loaded generator, the embedded magnets structure has also a much lower torque ripple. The permanent magnet is protected from pulsating fields by the iron of the rotor. Consequently, eddy currents are induced mainly on the surface of the magnet and in the region of the upper corners.

-0.9 0.0 0.9 Fig. 3: Magnitude of current density (A/mm 2 ) calculated using 3D model.

Lower losses in the magnets and a smaller cogging torque are possible by adding magnetic wedges between the stator teeth. The best effect is obtained using wedges made of stainless martensitic steel, which characteristics are similar to the magnetic material of the core. For composite wedges of low permeability the cogging torque is higher. Figure 3 shows the distribution of the axial component of the current density for the motor with wedges made from the stainless martensitic steel. Through the use of magnetic wedges made of martensitic steel it was possible to reduce the losses in the structure with the surface magnet to a value of 22.9 W, and for the structure with embedded magnets to 0.72 W (1.07 W for the structure without the wedges).

III. CONCLUSIONS

The paper proposes a methodology for computational assessment of the influence of the type of the material used for magnetic wedges and conductivity and segmentation of the permanent magnets on the power losses in the magnets. Additional analysis is provided and comparison of two designs, the surface mounted permanent magnets and the case of the magnets built into the rotor is presented. The comparison consider the torque ripple, which is important from the point of view of mechanical vibration, and power losses in the magnets, which may create a potential hazard due to local overheating of the magnets. A number of simulations, which enabled the dimensioning of the machine were carried out.

Fig. 1 :

 1 Fig. 1: Current density distribution on the surface of PM.

Fig. 2 :

 2 Fig. 2: Comparison of stator voltages for the embedded (1) surface magnet (2) (upper figure). Domparison of cogging torque for embedded (1) and surface magnets (-2). (lower figure).


Ruijun Liu

Zhang Bohua

Yu Zhuo

Han Zhiyuan

Wu Jingbang

Wang Ting

Vibration-based Terrain Classification Recognition Using a Six-axis Accelerometer

Keywords: vibration, terrain classification, Gabor atomic decomposition, time-frequency analysis

Aging is a key problem in our daily life which demands more technology such as the terrain classification recognition for the detailed path planning. In this article, a six-axis accelerometer was placed on the base of a walking stick to capture walking gait vibration information. The signal features in the collected information were used to reflect the terrain. We used a matching tracking method to decompose the signal data into linear combinations of Gabor atoms and a joint time-frequency analysis method to perform analysis on the topographic features in a three-dimensional space. Three types of terrains were classified from space-time data. Results show that the proposed algorithm could achieve an average accuracy rate of 88.6%, which is 4% higher than conventional methods.

1. INTRODUCTION

Nowadays, ageing problem has become increasingly serious with the changes in population structure and the significant increase in life expectancy [1]. In particular, elderly people who are ill and have mobility problems should be aware of the detailed terrain of the walking route before going out to avoid high-risk walking movements, such as uphill, bumpy or downhill terrain features. In order to achieve that, the detailed terrain data should be collected and classified beforehand. In this paper, we utilize a high reliable, small sized and facilitate integrated six-axis accelerometer, which is a MEMS-based sensor and can meet characteristics of autonomy and portability. The sensor is tied on a walking stick to collect the terrain data while people walking. After the data is collected, we extract different acceleration parameters to help identify these terrains. In order to collect the gait acceleration signal, we set up a device combining the accelerometer and the storage device, and input the signal into the algorithm model to extract the terrain category step by step, and finally get the experimental results. Specifically, we use Wigner-Ville Joint Time Frequency Analysis (JTFA) to study triaxial accelerometer signals (the time-domain and frequency-domain extended signals are often used to analyze nonstationary time series [START_REF] Ryu | Smart staff design to assist walking of elders[END_REF]), and to extract the distribution of JTFA features by principal component analysis (PCA) which is used to train and classify the effects of different terrains on gait behavior. Through the proposed method, the recognition accuracy and robustness of the algorithm can be improved.

At present, the methods of detection of terrain analysis can be mainly divided into two types. The first type of method uses human gait information or terrain distribution image obtained by video frames. In [1]- [START_REF] Chen | Multi-gait identification based on multilin-ear analysis and multi-target tracking[END_REF], different types of topographic relief features or motion track information features were extracted from surveillance video. These functions were used for terrain detection and analysis. The method of detecting terrain based on video frames has some difficult problems to overcome, such as the limitation of the number of cameras, installation positions, and camera capture views, which greatly limit the practical application of such method. The second type of method uses movement data collected and analyzed by wearable devices and related detection devices. [START_REF] Sun | Gait Characteristic Analysis and Identification Based on the iPhone's Accelerometer and Gyrometer[END_REF]- [START_REF] Rampp | Inertial sensor-based stride parameter calculation from gait sequences in geriatric patients[END_REF] utilized explore gravity acceleration sensors or gyroscopes for terrain detection. [START_REF] Sun | Gait Characteristic Analysis and Identification Based on the iPhone's Accelerometer and Gyrometer[END_REF] used iPhone-based accelerometer and gyroscope to analyze gait characteristics of different terrains. In [START_REF] Nandy | A study on gait entropy image analysis for clothing invariant human identification[END_REF], the subjects' gait data was acquired by a Wall-mounted motion recognition system, and machine learning algorithms were used to identify the gait to derive terrain features. Tong et al. [START_REF] Kozlowski | Deep neural networks for Terrain recognition task[END_REF] proposed a method based on Hidden Markov Model (HMM), which used a six-axis acceleration of the human body to detect terrain. David et al [START_REF] Lishani | Human gait recognition based on Haralick features[END_REF] introduced the design and implementation of hardware and software for a low-cost wearable system to monitor and retrieve information of human activity. In [START_REF] Singh | Vision-Based Gait Recognition: A Survey[END_REF], the authors used inertial sensors (3D accelerometer and gyroscope) carried by a trunk and thigh wearable system. The system operated on a generated countermeasure network (BAN), which makes the users more comfortable for long-term wear. Birdal [START_REF] Birdal | Human Identification Based on Gait Analysis: A survey[END_REF] reported that a key factor to analyze gait video is the combination of acceleration information in different dimensions. According to that, we can extract health features from 3d gait data to identify different human gait patterns. Some of the above methods require people with limited mobility to carry extra sensors to record their walking which increases their difficulty in walking. In addition, the related works in [START_REF] Zhao | Multisensor Image Fusion and Enhancement in Spectral Total Variation Domain[END_REF]- [START_REF] Lu | Deep Fuzzy Hashing Network for Efficient Image Retrieval[END_REF] involve Artificial Intelligence (AI), which are of great significance.

The study of gait analysis in this paper mainly focuses on the recognition of gait in different terrains.

The analysis is based on collecting walking data in a way that does not add additional burden for walking prediction. People with limited mobility usually uses walking sticks. Our portable sensors are installed on the walking sticks to record walking data. Next, based on the characteristics of this data, a machine-learning classifier is used to make behavioral predictions.

3. ABNORMAL BEHAVIOR CLASSIFICATION METHOD

Data collection

In the proposed algorithm, we connect a six-axis acceleration module to the Arduino UNO which is a microcontroller board, to measure and store the walking movement data (XYZ three-axis acceleration data). The data is recorded during a user's walking, and then transferred to the computer through a data transmission module. Next, a Kalman filter fusion algorithm is used to process acceleration data on the PC side. Finally, the collected data is integrated and classified using the Gabor atomic decomposition algorithm to determine the range of acceleration parameters corresponding to different terrains.

Data preprocessing

The embedded six-axis accelerometer described in Figure 1 can collect the acceleration data caused by human motion of X, Y and Z axes respectively. For different motor behaviors, such as uphill, downhill, and upright walking, the acceleration data is completely different.

To minimize the error caused by the initial installed position of the accelerometer, the amplitude of the signal vector is defined as formula (1):

𝑆𝑉𝑀(𝑡) = 𝑋 ! (𝑡) + 𝑌 ! (𝑡) + 𝑍 ! (𝑡) (1)
X(t), Y(t), Z(t)
represents the acceleration of the three axes respectively.

Before the SVM feature extraction, the masking threshold constraint Kalman filter is used to enhance the acceleration signal. The sample process of the n-th cleaning acceleration signal s(n) is simulated as follows:

𝑠(𝑛) = 𝑎 ! 𝑠(𝑛 -𝑖) + 𝑢(𝑛) ! !!! (2)
Here 𝑝 is the sample number,𝑎 ! is the i-th coefficient, 𝑢(𝑛) is the residual signal, and s(n) is assumed white noise with variance 𝜎 ! ! . As shown in Figure 2, after data preprocessing, the accuracy of the data is effectively enhanced, which is more in line with human gait characteristics. In particular, in the case of low SNR, the readability of the data is better.

Feature Extraction

During the training, it is necessary to collect and understand the vibration signal characteristics of the three types of terrain involved in this paper. Therefore, during the experimentation the user has to walk with the stick at six different speeds on three preset terrains to obtain the vibration information. The walking stick was equipped with a 100 Hz six-axis accelerometer. In order to facilitate the processing of experimental data, the collected data was divided into several segments. We can acquire 6 by 100 vectors per second, and the data was labeled according to the different terrain.

For the sake of data processing, we standardized the data, and normalized the vector in each direction of vibration to a vector with a mean of 0 in 100 dimensions and then took the Fourier transform.

Standardization reduces the influence of incorrect data on the final result, which is beneficial to the accuracy of the final result. After the normalization of the data, each indicator is in the same order of magnitude, which makes the evaluation results more comprehensive and improves the accuracy.

Common feature extraction algorithms, such as HOG (Histogram of Gradient) feature extraction algorithm [START_REF] Nandy | A study on gait entropy image analysis for clothing invariant human identification[END_REF] and SIFT (scale-invariant feature transform) local feature detection algorithm [START_REF] Rampp | Inertial sensor-based stride parameter calculation from gait sequences in geriatric patients[END_REF], have good performance in vibration feature extraction, but under the time efficiency, there is a large space for improvement. In this paper, a Gabor atomic feature extraction algorithm based on MP is introduced.

MP is a nonlinear iterative process, and the signal is decomposed into linear sums of known Gabor functions [START_REF] Rampp | Inertial sensor-based stride parameter calculation from gait sequences in geriatric patients[END_REF]. The decomposition of MP is realized through algorithm iteration. The number of iterations is determined by the residuals of the data. Firstly, the initial function value G selects (D) from the data to have the maximum function value. Next, the difference between the characteristic signal and the input signal is calculated. In a continuous iteration, the value of the Gabor function with the maximum product of the new given signal is selected. This process is repeated until the signal reaches target accuracy. The above steps can be defined as:

𝑅 ! 𝑥 = 𝑥 𝑅 ! 𝑥 = 𝑅 ! 𝑥, 𝑔 !" 𝑔 !" + 𝑅 !!! 𝑥 𝑔 !" = 𝑎𝑟𝑔𝑚𝑎𝑥 𝑔 !"∈! 𝑅 ! 𝑥, 𝑔 !" (3)
where 𝑥 is the signal atom. In the n-th alternations of MP, the atomic iteration signal 𝑥 is expressed as 𝑔 ! , and the difference of n-th is regarded as 𝑅 ! . The number of alternations is 𝑀. The final 𝑥 iteration result can be expressed by the following formula:

𝑥 ≈ 𝑅 ! 𝑥, 𝑔 !" 𝑔 !" !!! !!! (4)
Gabor function is used to analyze the time and frequency of different acceleration signals:

𝑔 !" = 𝐾 𝛾 𝑒 !! !!! ! 𝑐𝑜𝑠 2𝜋𝑓 𝑡 -𝑢 + 𝜑 (5)
𝛾 = (𝑢, 𝑓, 𝑠, 𝜑) is the time-frequency parameter. There are many factors that can affect it, including displacement, frequency, velocity and amplitude of vibration.

After summation of the WVD (Wigner-Ville Distribution) of the selected Gabor function, the time-frequency distribution of the acceleration signal can be obtained:

𝐸𝑥(𝑡, 𝑓) = 𝑅 ! 𝑥, 𝑔 !" ! 𝑊𝑔 !" (𝑡, 𝑓) ! !!! (6)
Then, the WVD matrix with time-frequency distribution is extracted as the feature by the algorithm in order to realize the recognition and classification of different terrains.

Currently, vibration-based terrain classification methods are only based on vertical (z-axis) detection.

The main reason is that the vertical acceleration is most sensitive to the terrain change and the change is the most obvious one. In order to further improve the recognition accuracy, we introduced the measurement of horizontal acceleration, that includes forward and backward (x-axis) direction, and left and right (y-axis) direction. By introducing these two dimensions of data, we can find the best feature recognition results in multi-dimensional data and different terrain.

In this paper, we collected the vibration data during a period of 1.5 minutes under three different terrains. According to the above method, 90 sets of data were collected for each terrain, thus 270 sets of data were collected in total. During the training process, 80% of the data were randomly selected for training, and the rest 20% of the data were used to test the accuracy of classification. The training process will be divided into two parts shown in Figure 3. In the training process of the model, we firstly preprocessed the collected data for noise reduction, and then used the algorithm mentioned above to obtain the gait features after fusion as input to the BP neural network for learning.

BP Artificial Neural Network Training

It is worth noting that because different walking sticks users may have different walking characteristics, the neural network can be used to effectively learn terrain features without being disturbed by the behavior patterns of different users. Later, in the test part, we used the walking data as the input system to get a three-way classification output as the terrain classification.

In order to classify three kinds of terrain, we also designed a three-layer BP neural network. The input layer is made up of 300 neurons in order to adapt to the corresponding input vector. The hidden layer contains 18 neurons, which corresponds to three different terrains the experiment, so there are three neurons in the output layer. After the training process, the neural network can carry out actual classification and collect unknown signals based on the experimental platform. The signal contains the 1×300 vector generated by the acceleration signal detected in three directions per second. After normalization and standardization, features were extracted with Gabor and then classified by trained neural network, and finally the evaluation results of terrain types were obtained.

4. EXPERIMENTS

Methods

The terrain classification experiment system includes a walking stick, a JY61 serial six-axis accelerometer/gyroscope, an Arduino board, a power supply (6v) and a switch, a USB transmission device, a 32G storage device, and a laptop computer. As shown in Figure 4, the total weight of walking stick and equipment is 0.76kg. A six-axis accelerometer/gyroscope is installed near the bottom of the walking stick. Acceleration data is automatically collected and stored during the user's walking. After that, the experimental data is read and imported into a computer for further analysis. The experimental terrain is divided into three types: straight cement pavement, inclined 30° uphill, and inclined 30° downhill. The schematic diagram of the experimental terrains is shown in Figure 5. We counted the walking speed using a speedometer in the experiment. The user was asked to walk back and forth uniformly on the three types of terrain at six speeds respectively, i.e., 0.6m/s, 0.8m/s, 1.0m/s, 1.2m /s, 1.4m/s, and 1.6m/s.

Results

When the walking stick moves on the three types of terrain, different acceleration signals will be generated dues to the different interactions between the walking stick and the ground. As shown in Table 1, we can see that the triaxial accelerometer data differs more on straight road.

The average value of the X-axis acceleration approaches 0, which is relatively gentle. The peak and range display show that the X-axis acceleration has a smaller vibration amplitude. The Y-axis acceleration has a negative mean value, and the highest peak value and range error, which indicate that the acceleration in the Y-axis direction changes greatly even when the road surface is gentle. The average value of the Z-axis acceleration is positive, the peak value and range are small, which indicate that the Z-axis acceleration changes smoothly on a smooth road. In the uphill slope of 30 degrees, the acceleration in the X-axis direction changes significantly, the average value is large and a significant peak appears and the peak in the Y-axis direction is positive, which indicate that the walking stick swings greatly in the Y-axis direction, and the acceleration in the Z-axis direction is still stable. In a downhill slope of 30 degrees, the average acceleration in the X-axis direction is negative, the peak value is low, and the range is enlarged; the acceleration fluctuation in the direction is huge, and it is obviously close to the same direction; the acceleration of the Z-axis acceleration is violent, the range is enlarged and the direction.is consistent. When downhill, it shows that the Z axis acceleration vibration amplitude is large and the direction is the same. It is not difficult to find through the chart.

Each type of terrain generates a unique acceleration signal, allowing the computer to identify each type of terrain and then process the data before classifying it with an algorithm. It can be seen from Figure 6 that the average accuracy rate of classification of a single sensor data is between 50% and 80%. Overall, the average accuracy rate of classification in the z-direction acceleration is slightly higher than that in the x-direction and y-direction. But it is always lower than the average accuracy of the classification of the three sensors. Because the fusion of 3 sensors can make use of the advantages of each sensor thereby improving the average accuracy of classification. Average classification accuracy(%)

0%

Walking speed(m/s) By calculating the average classification accuracy, we compared the algorithm with the SVM algorithm.

X-axis acceleration Y-axis acceleration Z-axis acceleration

Results after fusion

At medium and high speeds, the accuracy of our algorithm can reach more than 90%, which is higher than the SVM algorithm, as shown in Figure 7. At the same time, it can be seen that at the six speeds, the average accuracy of our algorithm finally reached 88.6%, which is about 4% higher than the traditional SVM algorithm.

5. CONCLUSION

In this paper, a terrain feature analysis method based on a six-axis acceleration sensor is designed and implemented. The detected acceleration signal is extracted and analyzed with Gabor atomic features, so as to classify different terrain in a three-dimensional space. The average classification accuracy of the experimental platform is 88.6% under six types of speeds. The experimental results prove that the accuracy of classification can be improved by analyzing the data acquired from the accelerometer. In the future work, will consider more sensors to obtain data, and apply the latest artificial intelligence methods to analyze the experimental results. Average classification accuracy(%)

6. ACKNOWLEDGMENTS

Walking speed(m/s)

After fusion SVM

Figure 1 .

 1 Figure 1. Terrain recognition device based on accelerometer

Figure 2 .

 2 Figure 2. Comparison of acceleration before and after data denoising

Figure 3 .

 3 Figure 3. Flow chart of terrain classification algorithm based on vibration

Figure 4 .

 4 Figure 4. Instrumented walking stick.

Figure 5 .

 5 Figure 5. The terrain kinds in experiments

Figure 6 .

 6 Figure 6. Comparison of classification result from different sensors' data at six walking speeds.

Figure 7 .

 7 Figure 7. Comparison of classification results.

Table 1 .

 1 Different terrain and corresponding gait data (1m/s)

	speed 1m/s		Straight road		Slope 30 degrees uphill	Slope 30 degrees downhill
	Direction of acceleration	mean value	peak	range	mean value	peak	range	mean value	peak	range
	X axes	0.07	0.63	1.26	0.37	01.23	1.46	-0.12	0.78	1.67
	Y axes	-0.12	-1.12	2.21	-0.26	-1.13	1.32	-0.56	-1.33	2.62
	Z axes	0.11	0.8	0.96	0.07	0.35	0.67	-0.23	-1.15	2.13

  This research was funded by the Beijing Natural Science Foundation grant number (4202016), National Natural Science Foundation of China grant number (62076012, 61972010), Research Foundation for

	Youth Scholars of Beijing Technology and Business University, PXM2020 014213 000017,
	QNJJ2020-27.

REFERENCES

[1] Berns, K. and Mehdi, S. A., "Use of an Autonomous Mobile Robot for Elderly Care", Advanced Technologies for Enhancing Quality of Life, pp. 121-126, 2010.