Turbulent shear flow of solids

The term"solid-state turbulence"may sound like an oxymoron, but in fact it is not. In this article we demonstrate that solid-state turbulence may emerge owing to a defining property of the solid state: the ability of a solid to retain its shape. We consider shear flow of layers of solids with different stiffness and show that the stiffer ones may spontaneously decompose into a set of blocks. This breakdown of isometry is a key to plasticity of solids and is fundamental for the occurrence of solid-state turbulence. To visualise the piecewise isometric transformations of the blocky structure in a turbulent flow regime, we use a heuristic model based on discretisation of a continuum into interacting"particles". The outcomes of the numerical experiments conducted support the occurrence of pulsations of velocity and pressure in plastically deforming solids and the emergence of vortices characteristic of classical turbulence. This phenomenon may have important practical implications for solid-state mixing as an ecologically beneficial alternative to conventional metallurgical processing routes.


Introduction
A flow which exhibits irregularities, chaotic pulsations of velocity and pressure, and random vortices of different scale that lead to mixing of the constituent substances is referred to as turbulent [1]. Almost any turbulent flow of a fluid is of shear character, its layers moving past each other under a shear force. Therefore, turbulence in shear flows is the main problem considered in that context [2]. Of special interest is the case of free shear flow, i.e., the flow of a fluid confined within another fluid.
It is almost a platitude to say that turbulence is one of the greatest enigmas of Nature, but this is, indeed, an apt statement [3]. Commonly, turbulence is associated with the gaseous or liquid aggregate states only. However, in metals deformed by severe plastic deformation (SPD) [4] a phenomenon clearly showing features of turbulence was observed. Similar effects are found in near-surface layers of metals engaged in dry friction [5,6]. In what follows, we shall refer to this phenomenon as solid-state turbulence (SST).
We would like to mention that as early as in 1953, Cottrell used the term 'turbulence' to describe a transition from the 'laminar' dislocation motion on a single slip system in stage I of strain hardening of 3 metals to multiple slip when several slip systems get activated in stage II [7]. Unlike this usage of the term as a figure of speech, we shall be talking about turbulence in a genuine sense, as a phenomenon characterised by chaotic changes in the flow velocity and the occurrence of unsteady vortices.
The article [8], in which a vortex structure in titanium processed by twist extrusion was documented, appears to be a first reported case of SST induced by SPD. The occurrence of the vortices was associated with rotational motion of the material under simple shear at the boundaries of the screwshaped part of the twist-extrusion die. The author of [9] invoked an analogy of SPD-induced SST with turbulence in fluids and hypothesised that vortices in a deforming solid are generated by obstacles hindering shear flow. Recently, several articles dealing with SST were published, but literature on the subject is still scarce. Thus, it was proposed to employ SST for mixing of metals aimed at fabricating new alloys or architectured materials. The principal SPD technique used to that end was high-pressure torsion (HPT) [10]. The process involves shear flows of two or more metals giving rise to vortices at their interfaces, which lead to random oscillations of the flow velocity and vigorous mixing of the metals, cf., e.g., [11][12][13]. The flow pattern and the shape of inclusions of the harder of the constituent metals within the relatively soft matrix resemble the picture observed in free shear flows of fluids. This was the motivation for computational modelling of the emergence of turbulence in a shear flow of multilayered metallic laminates in terms of the rheology of nonlinear viscous fluids [14] by finite element analysis [15,16] and molecular dynamics simulations [17,18]. Establishing a quantitative measure of the efficacy of mixing under SST, proposed and substantiated in [19], is also to be mentioned in this context. Although the possibility of the occurrence of solid-state turbulence may be evident from these reports, it is fair to say that the scientific community is largely unaware of the existence of the phenomenon. At any rate, the mechanisms of SST are not understood, and it appears timely to address them.
Investigations of SST as a stochastic nonlinear phenomenon is of great scientific significance, as it has certain universality and can occur at different length scales: from laboratory experiments to processes in the Earth's lithosphere [20]. The phenomenon is also of practical importance, as SST opens an avenue for solid-state mixing and physicochemical reactions under high pressure, whose potentialities in materials engineering can hardly be overestimated [13,21]. Indeed, SST offers an ecologically advantageous alternative to conventional metallurgical processing.
In this article, we propose and corroborate a hypothesis on the nature of SST and investigate some salient features of this phenomenon. To that end, a computational heuristic model was developed and applied. Fermi, Pasta, and Ulam appear to be the first to use a computer as a heuristic tool for solving a complex physical problem [22]. In the work that followed, Ulam and Pasta broadened the applications of this approach. Nowadays heuristic computer models are widely used to solve complex problems in a variety of research fields [23,24].
Working with a heuristic model is like playing an exciting game in which the researcher can observe the behaviour of a complex system in different situations created at will. This makes it possible to hypothesise on the origin of the phenomenon in question, identify the major governing parameters and their interplay, and get insights in the underlying physics, which enables development of quantitative models. For the 'game' to be efficient in this regard, the heuristic model must account for the salient features of the phenomenon meaning that it must capture the main governing factors. Accordingly, we start by establishing the factors that determine SST. 5

The major factors defining SST: modelling preliminaries
Since SST is found for solids of any chemical composition and size, we consider the principal defining characteristic of the solid state -the ability of a solid to retain its shape. Let us express this property in a formal way, so that it can be implemented in a model.
To be specific, we first look at crystalline solids. Their atomic lattice, which is responsible for long-range order in the system, is what distinguishes them from liquids. The existence of the crystal lattice owes to the physics of the interatomic interactions. We take this as a departure point of our analysis and pose the problem of deformation of a solid in geometrical terms. This will allow us to draw some rather general conclusions regardless of the chemical or phase composition of the solid.
We denote the initial and the final configurations of a solid body by A and B, respectively. The geometrical operation which transforms A to B is denoted by G [25]. It is known that plastic deformation has no substantial effect on the lattice parameter, changing it by no more than ~0.1%.
This means that if G is to represent plastic deformation, it can be considered to belong to the class of isometric transformations. These are transformations under which the distance between any two points in the body remains unchanged [26]. Transformations of this kind include rigid displacements, rotations, and mirror reflections, see Fig. 1. If the solid does not undergo plastic deformation, the transformation G is the same for its entire volume, so that the configurations A and B differ only in the spatial location and/or orientation. In the opposite case, G must include discontinuities of displacements, rotations, or mirror reflections, which puts it in the category of isometric transformations with singularities, commonly referred to as piecewise isometries [27].
So far, we considered only crystalline materials, but the conclusion that the geometry of a solid undergoing plastic deformation changes by means of piecewise isometries has general validity and applies to amorphous solids, as well. It also applies to fracture processes, which occur by subdivision of a solid into fragments that are isometric to their original regions within the solid before fracture.
Comminution of a rock is an example of such a process. 7 The concept of piecewise isometry has been applied to describe various phenomena and systems, such as, e.g., pattern formation in Nature [28], mixing of powders [29][30][31], and behaviour of dynamic systems [27]. The concept of plastic deformation and fracture of solids being describable in terms of piecewise isometries provides explanations for a whole range of observed effects. As an example, isometry discontinuities associated with the deformation of solids have a stochastic character and can give rise to mass transport and efficient mixing of co-deforming solids. We believe that isometric discontinuities may also be one of the possible reasons for anomalously fast diffusion in severely deformed metals [37].
By the same token, we see the origin of SST in piecewise isometries accompanying deformation and ruptures of a solid. The maximum distance between singularities of isometry determines the spatial scale of SST. The characteristic size of the isometric fragments of a solid is governed by the constraints put on the flux of the material. For example, if a solid is sheared between two rigid surfaces, it will imitate laminar flow, the deviations from laminarity becoming more and more evident as the length scale of the observation is decreased. If the constraining boundaries are not rigid, the non-laminar flow will disturb them more and more as the amount of shear increases, so that turbulence will become manifest at macro scale.
We believe that it is a tight constraint on deformation in HPT that explains why it results in the smallest grain size (i.e., the greatest degree of fragmentation) among the SPD techniques [4].
Thus, the first and foremost tenet of a prospective heuristic model of solid-state turbulence is the assertion that the flow geometry of a solid body is associated with piecewise isometries. A second important ingredient of the model is the establishment of conditions leading to shape variation of a solid, i.e., to the generation and propagation of isometry discontinuities. In solid mechanics, such conditions correspond to plasticity and fracture criteria. A rather general criterion of this kind is that of a critical magnitude of the elastic strain energy, which is related to the yield stress under shear deformation.
In the next section, we present our mathematical model of solid-state turbulence. The model may be somewhat naïve compared to traditional models of plasticity and fracture, yet it will be shown to have the advantage of accounting for the salient features of SST.

The model of solid flow
Our aim was to employ a model of plastic flow in solids that would be robust and simple, yet capable of capturing the above principal factors defining SST. Based on such a model, we endeavoured to demonstrate through numerical experiments that it is the patchiness of the emerging geometrical patterns that breaks the laminarity thus giving rise to turbulent flow. The patchiness, in turn, is contingent on the degree of order in the system. We are obviously dealing here with a particular case of a general problem: how do meso-and macroscopic properties (including the elastic and plastic ones) 9 and the associated geometrical patterns derive from dynamic processes described at a deeper theoretical level? Suitable modelling tools are provided by discretisation of a continuum and tracing the dynamics of 'particles' representing it. For example, the model proposed in [38] in terms of the dynamics of discrete 'particles' yields a broad range of mass and energy transport regimes -from the solid to the gaseous state -for the simplest case of a repulsive interaction potential between the 'particles'. This approach provides control of the correlation radius between the elements of an effective continuum, i.e., the degree of order, which is crucial for the occurrence of SST. To avoid confusion, it should be stated from the outset that the use of the notion of 'particles' here and below does not imply that individual atoms are considered. Recently, а significant progress with understanding structural lubricity in soft and hard matter systems was achieved by using a similar modelling approach [39]. At present, it is still not possible to advance to the satisfactory length and time scales by using classical molecular dynamics that operates with atoms or molecules. That is why modelling approaches in which representative larger-scale objects whose motion in its entirety reflects the pertinent aspects of the behaviour of the system have come to the fore. These include, for example, the method of moveable cellular automata [40], smoothed particle hydrodynamics [41], and the discrete automata method [42,43]. Not only do these approaches show a qualitative accord between modelling and experiment, but they also enable a good quantitative agreement for various systems considered. However, they are computationally costly. Besides, visualisation and processing of the results of the computations are very involved [44].
Our concern in the present work was primarily the extraction of qualitative information from the numerical simulations, with a convenient access to direct visualisation for various scenarios and over large periods of time. As a method of choice, we used the isotropic 'moveable automata' technique (cf., e.g. [45]), which involves a weak long-range attractive interaction between particles in addition to a short-range repulsive one. The form of the interparticle interactions used ensures the existence of a minimum of the overall potential and the attendant equilibrium state. Such a system does not require 'walls' created by the boundary conditions to form compact fragments of a crystalline lattice. Even for the case when the average density of the material is lower than that required for contiguous space filling, the particles show a tendency to aggregation. The freedom of the solid to deform plastically is still retained, as the lattice is subdivided in fragments capable of moving over distances far greater than the lattice constant. In doing so they can rotate and deform elastically nearly as a whole, or fuse with other fragments that may have been distant initially. The ability to rupture enables the lattice to engulf particles of other substances that were separated from it. Naturally, the combinations of attraction and repulsion are different for the sub-systems of a two-component or, generally, multicomponent system, as well as between the subsystems. It can be shown that a strong mutual repulsion promotes separation of subsystems, whereas an attraction enhances the tendency to mixing.
We used this observation in connection with modelling of the behaviour of materials with different stiffness that show some (but not exceedingly strong) proclivity for separation. The details of the model used are given in the Supplementary Information.
In solving the equations of motion, Supplementary Information Eq. (S3), we confined ourselves to problems where the thermal energy , where T is the absolute temperature and kB is the Boltzmann constant, is negligibly small compared to other energy terms involved. It can be proven [38], that in the limit of high kinetic energy density, Ekin >> ij C , where ij C is a stiffness parameter associated with the magnitude of the attractive interaction between particles i and j, the system described by Supplementary Information Eqs. (S1) and (S2) behaves like a gas of nearly free-flying particles. In the 11 opposite limit case, Ekin<< ij C , a crystal lattice is formed spontaneously. The lattice constant a is determined by the equilibrium distance between the particles. The nodes of the lattice oscillate around the bottoms of the potential valleys formed by the neighbouring particles. This oscillatory motion can be described in terms of a harmonic Hamiltonian slightly perturbed by non-linear terms. In this sense, the system deforms elastically. While in the gaseous state the system is isotropic on average, in the solid state its lattice has a hexagonal symmetry if the interaction potential between the particles is isotropic. In both limit cases, the system exhibits a regular dynamic behaviour. However, due to a difference in symmetry, there is a transition between the two states via a mixed disordered state, which, in the presence of dissipation, represents a viscous fluid. Thus, the minimalist model considered accounts for the occurrence of all three aggregate states of the system.  It should be noted that such geometrical constructions are helpful both for qualitative presentation and for quantitative evaluation of the numerical simulation results, e.g., in the form of uniquely determined histograms for particle distributions.

Solid-state turbulence as revealed by computational experiments
Building upon the knowledge of the basic properties of the model discussed above, we can now turn to considering turbulence in a composite system comprised by materials with different stiffness. We construct the initial system as a layer of a stiff material sandwiched between two layers of a more compliant material. The latter are set in motion by two rigid plates shifted horizontally in the opposite directions with the same constant velocity, . V const  Figure 3 shows the results of numerical simulations. Movie_02.avi, clearly illustrates that at this initial stage the movement of the external plates is transmitted only to the more compliant sub-system, the stiff one being practically indifferent to this movement. Since the sandwich structure was created naturally, through spontaneous ordering of the layers of two different kinds before shear deformation, the interfaces between the layers exhibit some roughness. A phase separation between the layers is furnished by mutual repulsion between the particles of the two kinds considered to be the same as that between the particles of the same kind, combined with a smaller (or even non-existent) attraction. The smoothness of an interface and the separation distance between the layers of different kinds (or the opposite -the interdiffusion of the two particle species) are controlled by the magnitude of the attractive interaction.
To be specific, here we confine ourselves to an intermediate case when a tendency to phase separation does exist but is not excessive. Such a case is shown in Fig. 3c. Solid-state turbulence promotes mixing of the two materials of which the sandwich structure is composed. The degree of mixing can be quantified using a criterion proposed in [19]. It is based on the use of the mixing index The minimum value, 2 min GC  , is attained for a uniform distribution, C = const, while a maximum is reached when the concentration at each node of the grid is equal to zero or unity. Hence, the index Q can be used to quantify the degree of mixing of the two components. We used this procedure to characterise the disruption of a stiff layer confined by two compliant layers described above. Figure 4 shows a snapshot of the system taken at an intermediate stage of its shear flow, along with the time 17 dependence of the mixing index Q. It is seen that, as expected, the mixing index gradually drops with the progress of the shear flow process.

Discussion
The results of the numerical experiments conducted to model SST allow us to draw some general conclusions and establish the salient properties of free shear turbulence in solids.
As shown in Sections 3 and 4, long-range order in the system gives rise to a blocky structure of both materials of a 'sandwich'. In the shear flow, the blocks adjust to each other, get fragmented or agglomerate, but most of their bulk undergoes isometric transformations, rather than deforming. Hence, the model underlying the numerical simulations corroborates our hypothesis of piecewise isometric character of the deformation of solids. The success of the model with adequately emulating SST gives reasons to believe that this characteristic feature of the deformation of solids is the root cause of solid-state turbulence.
What features of a real experiment correspond to the fragmentation of solid flow into blocks? It can be said with certainty that these blocks are not microstructural entities of the material (such as grains or sub-grains). Indeed, the solid considered in the model is structureless. The blocks represent chunks of material that are carried by the overall flow without being deformed. Their characteristic size, and the spatial pattern of SST in general, is dictated by the length scale of the order in the system and reflects the level of freedom of motion of the material points. Computer experiments presented in Section 3 showed that an increase of the block size is associated with the growth of the yield strength of the model material.
When the shear flow of a solid is constrained by rigid boundaries, segmentation of the blocks down to sizes much smaller than the width of the shear layer occurs. Viewed at macro scale, the solid-state flow appears to be laminar, but at smaller length scales disruption of the laminar flow becomes evident. If the shear flow of a solid layer occurs in a relatively compliant environment, the big blocks into which it gets subdivided are less prone to further fragmentation than in the case when the layer is confined by rigid walls. In the latter case, the turbulent character of the flow is manifested at a length scale prescribed by the thickness of the layer, and SST is pronounced at macro scale. This macroscopic turbulent flow leads to intensive mixing of the rigid and the compliant constituents of the sandwich.
The deformation of a model material is caused by an external shear force, which is transferred into the bulk of the flux through interaction between particles. Since particles are agglomerated in blocks, an entire block acted upon by a force from its environment can be considered. According to the laws of mechanics, this force produces a central force and a momentum, which leads to block rotation. The model employed here accounts for this effect, which enables a qualitative description of the phenomenon of SST. To treat SST in a continuum mechanics model, gradient theories of plasticity would need to be invoked and a complex non-uniform stress-strain state of a solid under a shear flow would have to be considered [46]. 19

Conclusion
In this article we hypothesised that the phenomenon of solid-state turbulence is a result of the fundamental property of the solid state: the ability of a solid to retain its shape and to change it only under sufficiently high stresses. To validate this hypothesis, we employed a discretised, multi-particle It is opportune to represent the interaction potential by a pair of the Gauss potentials: The simplicity of the equations of motion, is deceptive. Summation at every time step over all possible neighbours, including very distant ones and the fictitious image particles stemming from the periodic boundary conditions, is required. In addition, in the case of contact with a thermal bath the boundary conditions need to be replaced with more realistic ones. If the temperature T of the thermal bath at a given boundary (or within the system) is non-zero, an array of random  -correlated Langevin forces need to be included. These are to satisfy the fluctuation-dissipation theorem in the form of ( , , and dt is the discrete time step. Interacting particles exchange momentum, and hence a dissipation channel acting to equilibrate the relative velocities of particles that happen to be within the distance v c close to the equilibrium one needs to be introduced. This is done by including an additive force on the i-th particle, with yet another dissipation constant . As a result, the equations of motion assume the following form: They can be integrated by using Verlet's method, which conserves the energy of the system at each time step, provided no energy is supplied externally through mechanical work or temperature variation. 27

Supplementary Information 2. Calculation of the force
Here we calculate a measurable quantity: the external force acting on the system. First, we calculate the work produced by the forces acting on each particle of the system or the corresponding power  From this dependence the yield force, or yield stress, can be identified. In our computer experiments, 29 this dependence resembled that for the resistance force for sliding a solid on a surface with random roughness in the stick-slip regime. This is not surprising, as in our case at any given moment at least two layers of the material separated by a rough interface are engaged in a relative frictional motion. In  It should be mentioned that the position of the maximum of a histogram, i.e., the most probable magnitude of the force, is practically stiffness-independent, while due to the tail of a distribution the average force is shifted to larger values when the stiffness is increased. As this shift owes to relatively seldom, but large, maxima of the force, one can conclude that the measured force is governed by sharp discontinuous stick-slip events, which are more pronounced in stiff systems. Figure S3. Results of computations of the time-averaged force as a function of the stiffness of the system (normalised with respect to a certain value 11 * , which was fixed in all numerical experiments).
The dash-dotted horizontal line corresponds to a background common to all systems considered.
It can be shown that the common background stems from the viscous contribution to the resistance to plastic flow, which in the model considered is present for all systems -whether stiff or soft ('fluid') systems alike. To exclude it, we conducted similar numerical experiments on soft systems, which were otherwise similar to those with a much lower dissipation constant of drop of the background level, which is to be subtracted from the force readings when determining the plastic yield limit, is determined for a fixed value of  , was observed.
In principle, the integral force can be determined both for a large system or, alternatively, through averaging over a large number of realisations. However, analysis of the time dependence of the force, even when averaged over a relatively small number of realisations, shows that it becomes rather flat.
For a prolonged steady state, the results of averaging over the histograms in Fig. S2 and the time average converge. This is in keeping with the 'ergodic hypothesis' which states that for a sufficiently large system the average over time and that over the realisations coincide. We used this observation to determine the average force for a representative large number of cases and to establish its dependence on the key parameter, 11 C , which governs the stiffness of the system. The calculated values of the force, which were shown to be in accord with the expected values, are compiled in Fig.   S3. The horizontal dash-dotted line corresponding to the background, which is to be subtracted from the data in the case of . const  

Annotations to the videos Supplementary Movie_01.avi
Typical motion of a solid under shear deformation. The colour code is the same as in Fig. 3. It is seen that initially the solid deforms elastically, the particles' motion smoothly following the displacement of the boundary plates. The corresponding force vs. displacement curve is smooth. When a critical shear is reached, the bonds between neighbouring particles start breaking. However, locally the lattice domains into which the solid gets fragmented move as a whole, switching to a different