Table Of Content

Is dislocation flow turbulent in deformed crystals? Woosong Choi∗ and Yong S. Chen, Stefanos Papanikolaou, James P. Sethna† Laboratory of Atomic and Solid State Physics (LASSP), Clark Hall, Cornell University, Ithaca, New York 14853-2501, USA Intriguing analogies were found between a model of plastic deformation in crystals and turbulence in fluids. A study of this model provides remarkable explanations of known experiments and predicts fractal dislocation pattern formation. Further, the challenges encountered resemble those in turbulence,which is exemplified in a comparison with theRayleigh-Taylor instability. 2 1 0 From horseshoes and knives to bridges and aircrafts, thatignoresmanyimportantfeaturesofplasticdeforma- 2 mankindhasspentfivemillenniastudyinghowthestruc- tion of crystals, including yield stress, work hardening, n tural properties of metals depend not only on their con- dislocation entanglement, and dependence on material a stituents, but also how the atoms are arranged and re- properties [1]. We focus on the complex cellular struc- J arrangedas metals are cast, hammered, rolled, and bent turesthat develop in deformed crystals,which appear to 8 1 into place. A key part of the physics of this plastic dis- be fractal in some experiments [2]. These fractal struc- tortion is played by the motion of intrinsic line defects tures are reproduced by our continuum dislocation dy- ] called dislocations, and how they move and rearrange to namics (CDD) [3] theory (see Figure 1a). h p allow the crystal to change shape. Not only do the resulting patterns match the experi- - Here, we describe the intriguing analogies we found mentalones,but the theoryalso has richdynamics, akin p between our model of plastic deformation in crystals to turbulence. This raises a question: Is the dislocation m and turbulence in fluids. Studying this model led flowturbulent? Here,wefocusonexploringthisquestion o us to remarkable explanations of existing experiments by building analogies to an explicit turbulence example: c and let us predict fractal dislocation pattern formation. the Rayleigh-Taylorinstability. As we describe, our the- . s The challenges we encountered resemble those in turbu- ory displays similar conceptual and computational chal- c i lence, which we describe here with a comparison to the lengesasdoesthisexample,whichreassuresusthatwe’re s y Rayleigh-Taylorinstability. on firm ground. h For brevity, we offer a minimal problem description This CDD model [3–5] provides a deterministic expla- p nation for the emergence of fractal wall patterns [3, 4] [ in mesoscale plasticity. The crystal’s state is described 3 by the deformation-mediating dislocation density ̺ij – v wherei denotes the directionofthe dislocationlines and 1 j their Burgers vectors [1] – and our dynamical evolu- 5 tion moves this density with a local velocity V , yielding 3 ℓ a partial differential equation (PDE): 5 . 5 4 ∂ ̺ −ε ∂ (ε V ̺ )=ν∂ ̺ (1) 0 t ij imn m nℓk ℓ kj ij 1 HereV isproportionaltothenetforceonit(overdamped 1 ℓ : motion), coming from the other dislocations and the ex- v ternal stress. That is, i X D r V = σ ǫ ̺ a ℓ |̺| mn ℓmk kn (a)ContinuumDislocationDynamics (b)Turbulence where σ is the localstress tensor,the sum of anexternal FIG. 1: Comparison of our continuum ext stress σ and the long-range interactions between dis- dislocation dynamics (CDD) with turbulence. locationisj. σint = K (r−r′)̺ (r′)dr′, with K (a) Dislocation density profile as it evolves from a ij ijmn mn ijmn the function repreRsenting the stress at r generated by ̺ smooth random initial condition. The structures form at r′ [1]. The term proportional to ν is the regularizing fractal cell wall patterns. Dark regions represent high quartic diffusion term for the dislocation density (an ar- dislocation density. (b) Rayleigh-Taylorinstability at a tificial viscosity), which we’ll focus on here. (In fact, the late time. The fluid (air) with two layers of different equation we simulate here is further complicated to con- densities mix under the effect of gravity. The emerging strain the motion of the dislocations to the glide plane flows exhibit complex swirling turbulent patterns. The whileminimizingtheelasticenergy[3,4].) Thedetailsof color represents density (red for high, blue for low). our equations aren’t crucial: dislocations move around 2 with velocity V~, pushed by external loads and internal 100 stresses to lower their energies. Our equation is nonlinear, and it’s exactly this non-linearity that makes our 10 theorydifferent frommoretraditionaltheories ofcontin- 2 uum plasticity. 1/> Turbulence is an emergent chaotic flow, typically de- Ν2|| scribed by the evolution of the Navier-Stokes equations ρ − N=16 at high Reynolds numbers: Ν N=32 2 ρ(∂t~v+~v·∇~v) = µ∇2~v+f~ (2) ρ<|| NN==61428 ∂ ρ+∇·(ρ~v) = 0 N=256 t 0.001 N=512 where ρ is the local density of the fluid with velocity ~v under the application of local external force density f~. 0.0001 The term proportional to µ is the fluid viscosity, and µ 0.01 Time 1 is inversely proportional to the Reynolds number. (a)Non-convergence indislocationdynamics(seeEq.(1)) Despite this Navier-Stokes equation’s enormous suc- cess in describing various experiments, there are many mathematical and numerical open questions associated 0.3 with its behavior as µ → 0. In this regime, complex scale-invariantpatterns of eddies and swirls develop in a 1/2 > way that isn’t fully understood: turbulence remains one 2 of the classic unsolved problems of science. Ν|| ρ Howisνrelatedtoµ? Eq.(2)canbewrittendifferently − N=16 Ν N=32 by dividing the whole equation by ρ, in this equation 2ρ N=64 µ/ρ (in the incompressible case) is called the kinematic <|| N=128 N=256 viscosity (usually denoted as ν). Our artificial viscosity ν in Eq.(1) is analogous to this kinematic viscosity. For 0.01 turbulence, µ in Eq.(2) is given by nature. In contrast, ourν inEq.(1)isaddedfornumericalstability;itsmears 0 1 2 Time 4 5 singularwallsto giveregularizedsolutions. This isphys- (b)Non-convergence inturbulence icallyjustifiedbecausetheatomiclatticealwaysprovides dynamics(Rayleigh-Taylorinstability) a cutoff scale. How do we know that this artificial term FIG. 2: Non-convergence exhibited in both gives the ‘correct’answer (given that there can be many plasticity and turbulence. As time progresses,the different solutions to the same PDE)? Numerical meth- curves, which are initially monotonically decreasing, flip ods for shock-admitting PDEs are validated by showing order and become non-convergent(where the lines cross that the vanishing grid spacing limit h → 0 gives the each other). Red arrows show where the convergence is same solution as the ν → 0 limit (that is, the viscosity lost. solution). Wewillarguethatbothourmodelofplasticity and the Navier-Stokes equations do not have convergent solutionsasµorν →0. We’rereassuredfromtheturbu- lenceanalogyandaresatisfiedwithextractingphysically version of Navier-Stokes) faster than its speed of sound. sensibleresultsfromourplasticitytheory–viewingitnot The sonic boom is a sharp jump in density and pres- as the theory of plasticity, but as an acceptable theory. sure, which causes the continuum equations to become To solve Eq.(1), we implemented a second order cen- ill-defined. Our PDEs, depending on gradients of ̺, be- tral upwind method [6] especially developed and tested come ill-defined when ̺ develops an infinite gradient at for conservation laws, such as Eqs.(1) and (2). The a dislocation density jump. method uses a generalized approximate Riemann solver The numerical methods we use are designed to appro- which doesn’t demand the full knowledge of characteris- priately solve the so-called Riemann problem: the evo- tics [6]. For the simulations of Navier-Stokes dynamics lution of a simple initial condition with a single step in (Eq.(2)) we use PLUTO [7], a software package built to theconservedphysicalquantities. Forhyperbolicconser- run hydrodynamics and magnetohydrodynamics simula- vationlaws,exactsolutionsof the Riemann problemcan tions, using the Roe approximate Riemann solver. be obtained by decomposing the step into characteristic Accurately capturing singular flows is a challenge in waves. However,inmostnonlinearproblems,findingex- computationalfluiddynamics. Aclassicexampleofsuch act Riemann solutions involves iterative processes th at singularitiesisthesonicboomthathappenswhenanob- are either slow or (practically) impossible, and thus ap- ject passes through a compressible fluid (described by a proximatesolutionsareemployedinstead. Bothmethods 3 we use are approximate in different ways, but are quali- For many problems, researchers have shown that tatively similar. adding an artificialviscosityand taking the limit to zero These sophisticated methods are designed to handle yields a weak solution to the problem. For some prob- the kind of density jumps seen in sonic booms. In our lems, the weak solution is unique, while for others there dislocationdynamics,though,wehaveamoreseveresin- can be several: different numerical methods or types of gularity that forms – a sharp wall of dislocations that regularizingviscositiescanyielddifferentdynamicsofthe becomes a δ-function singularity in the dislocation den- singularities. Thismakesphysicalsense: ifasingularde- sity (as ν →0). These δ-shocks are naturally present in fect (a dislocation or a grain boundary) is defined on crystals – they describe, for example, the grain bound- an atomic scale, shouldn’t the details of how the atoms aries found in polycrystalline metals, which (in the con- move (ignored in the continuum theory) be important tinuumlimit)formsharpwallsofdislocationsseparating for the defect’s motion? In the particular case of sonic dislocation-free crystallites. Unfortunately, the mathe- booms, the microscopic physics picks out the viscosity maticalandcomputationalunderstandingofPDEsform- solution (given by an appropriate µ → 0 limit), lead- ing δ-shocks is relatively primitive; there are only a few ingmathematicianstolargelyignorethe questionofhow analytic and numerical studies in one dimension (for an micro-scale physics determines the singularities’ motion. example, see [8]). Currently, to our knowledge, there’s However,ourproblemsherearemoreseverethanpick- no numerical method especially designed for δ-shock so- ing out a particular weak solution. Neither our disloca- lutions. Moreover,inastrictmathematicalsense,several tiondynamicsnorthe Navier-Stokesequation(withvery properties of Eq.(1) – nonlocality and mixed hyperbolic high Reynolds number) converge in the continuum limit andparabolicfeatures–haven’tbeenproventopermita even for gross features, whether we take the grid size to successfulapplicationofshock-resolvingnumericalmeth- zeroin the upwind schemes or we take ν →0 (or µ→0) ods. as a mathematical limit. In the simulations presented here, we won’t add an Figure2ashowsaquantitativemeasureoftheoursim- explicit viscosity (so µ = ν = 0); instead, we have an ulation’s convergence as a function of time, as the grid effective numerical dissipation [6] that depends on the spacing h = 1/N becomes smaller. We measure conver- grid spacing h as hn, where n depends on the numerical gence using the L2 norm method used. (Eq.(2) with µ = 0 is the compressible 1/2 Euler equation. Although we present our simulations as ||̺2N −̺N||2 ≡ k̺2N −̺Nk2dx small numerical viscosity limits of Navier-Stokes, they (cid:18)Z (cid:19) could be viewed as particular approximate solutions of d d these Euler equations.) where ̺2N has been suitably smeared to the resolution Figure 1 shows typical emerging structures in simu- of ̺N. (Normally we’d check the difference between the d lations of both the CDD Eq.(1) and the Navier-Stokes current solution and the true answer, but here we don’t dynamics Eq.(2): both are complex, displaying struc- knowthe trueanswer.) Here,We study the relaxationof tures at many different length scales; sharp, irregular asmoothbutrandomlychoseninitialcondition–thatis, walls in the CDD and vortices in Navier-Stokes. Al- a perfect single crystal beaten with mesoscale hammers though it might not be surprising to professionals in with round heads – as a function of time. We see that fluid mechanics that nonlinear PDEs have complex, self- for short times these distances converge rapidly to zero, similar solutions, it’s quite startling to those studying implying convergence of our solution in the L2 norm. plasticity that their theories can contain such complex- However,ataroundt=0.02to0.2,thesolutionsbeginto ity(eventhoughthiscomplexityhasbeenobservedinex- become increasingly different as the grid spacing h→0. periments [2]): traditional plasticity simulations do not This worried us at the beginning because it suggests lead to such structures. that the numerical results might be dependent somehow These richandexotic solutionsdemandscrutiny. How on the artificial finite-difference grid we use to discretize do we confirm the validity of our solutions? For contin- the problem, and therefore might not reflect the correct uum PDEs solved on a grid, an important problem that continuum physical solution. We checked this by adding needs to be addressed is the effect of the imposed grid. the aforementioned artificial viscosity ν in Eq.(1). We Traditionally,it’sexpectedthatasthegridbecomesfiner, foundthat it convergesnicely when ν is fixed asthe grid thesolutionislikelytobeclosertotherealcontinuumso- spacing goes to zero. However, this converged solution lution. Fordifferentialequationsthatgeneratesingulari- is not unique: it keeps changing as ν → 0. So, it’s our ties, one cannot expect simple convergence at the singu- fundamentalequationofmotion(Eq.(1))andnotournu- lar point! How do we define convergencewhen singulari- merical method that fails to have a continuum solution. ties are expected? For ordinarydensity-jump shockslike This would seem even more worrisome: How do we un- sonicbooms,mathematicianshavedefinedtheconceptof derstand a continuum theory whose predictions seem to a weak solution: it’s a solution to the integrated version depend on the smallest studied length scale (the atomic of the original equation, bypassing singular derivatives. size)? 4 (a)h=1/128 (b)h=1/256 (c)h=1/512 FIG. 3: Continuum dislocation dynamics. Simulation results at t=1.0 of our CDD Eq.(1) at different grid sizes (h), starting from a smooth initial condition. We use periodic boundaries in both horizontal and vertical directions, and all physical quantities are constant along the perpendicular direction. itydecreases(forfixedinitialconditionsandloading)not only do the small-scale eddies get smaller, but also the position of the large-scaleeddies at fixed time change as the viscosity or grid size is reduced. Figure 2b depicts the convergence behavior of a simulation of the Rayleigh-Taylor instability (as in Fig. 2a). The instability triggers turbulent flow, and like our dislocation simulations, convergence is lost after t ∼ 2.0 as the grid spacing h gets smaller. Our choice of the Rayleigh-Taylor instability for comparison is motivated bythepresenceofrobustself-similarfeatures(suchasthe “bubbles” and “spikes” in Fig. 4 [9]), and by the spatio- temporally non-converging features of the initially well- (a)h=1/128 (b)h=1/256 (c)h=1/512 defined interface. The interface between two fluid densities is analogous FIG. 4: The Rayleigh-Taylor instability of the to our dislocation cell walls. Even though the Rayleigh- Navier-Stokes equation. The Rayleigh-Taylor Taylor instability is different from homogeneous turbu- instability is a fluid mixing phenomenon that occurs lence in important ways, we also verified that the lat- when an interface between two different fluid densities ter shows similar spatio-temporal non-convergence but is pulled by gravity. These simulation results here are statisticalconvergence(simulating the Kelvin-Helmholtz the Rayleigh-Taylorinstability at t=4.0, for µ→0, at instability for compressible flow in 2D). The Rayleigh- different grid spacings (h) with periodic boundaries in Taylor instability provides the best visualization of the the horizontal direction and fixed boundaries along the analogy between the two phenomena, but this non- vertical. The initial condition has density interface with convergence appears to be more general. Fig. 4 shows a single mode perturbation in the vertical velocity. The density profiles at intermediate times for the Rayleigh- system size is (L ,L )=(1.0,2.0). x y Taylor instability. The figure shows the formation of vortices(swirlingpatterns),andthesimulationslooksig- nificantly different as the grid spacing decreases. Again, this is analogousto the correspondingsimulationsof our It’s here that the analogy to turbulence has been cru- dislocation dynamics shown in Fig. 3, where larger cells cial for understanding the physics. It’s certainly not ob- continuetodistortandshiftasthegridspacingdecreases. viousthatthelimitofstrongturbulenceµ→0inNavier- Stokes (Eq.(2)) should converge to a limiting flow. Ac- Ifthesimulationsaren’tconvergent,howcanwedecide tually, ourshortexperiencesuggeststhat there is no vis- if the theory is physically relevant and can be trusted cosity solutionfor Eq.(2). Turbulence has a hierarchyof to interpret experiments? In turbulence, it has long eddies and swirls on all length scales, and as the viscos- been known that, as vortices develop, self-similar pat- 5 theexistenceoffinite-timesingularitiesisatopicofactive research: eventhoughlocal-in-timeanalyticsolutionsare easily shown to exist, global-in-time analytic solutions can be proven to exist only for special cases, such as in the 2D incompressible genuine Euler equation [12]. In 3D,the mechanismofvortexstretchingis conjecturedto lead to finite-time singularities [13], even though there are still crucial open questions. Despite its complexity, turbulencecanbeconcretelystudiedinspecialcases. For our example of the Rayleigh-Taylor instability, the two- fluidinterfacegetsdistortedand“bubbles” form(Fig.4); overtime,thebubblesexhibitemergentself-similarchar- acteristics [9], showing statistical convergence. However, there is no spatio-temporal convergence, because the interface develops complex, turbulent features as the grid becomes finer (see Figs. 2b and 7). FIG. 5: Statistical properties and convergence Although the continuum dislocation dynamics (CDD) Sometimes science seems to be fragmented, with in- simulations with different grid sizes are dependent fields whose vocabularies, toolkits, and even non-convergent(Fig. 2a), the statistical properties are philosophiesalmostcompletelyseparate. Butmanyvalu- the same. The dislocation density correlationfunction ableinsightsandadvancesarisewhenideasfromonefield is plotted here for different simulation sizes at the same are linked to another. Computational science is provid- time, exhibiting consistent power laws. ing a new source of these links, by tying together fields that can fruitfully share numerical methods. Our use of well-established numerical methods from terns arise in the flow and exhibit power laws in the the fluids community made it both natural and easy to energy spectrum and in the velocities’ correlation func- utilize their analytical methods for judging the validity tions[10]. Asuccessfulsimulationoffully developedtur- of our simulations and interpreting their results. bulence isn’t judged by whether the flow duplicates an We thank Alexander Vladimirsky, Randall LeVeque, exact solution of Navier-Stokes! Turbulence simulations Chi-Wang Shu, and Eric Siggia for helpful and inspir- study these power laws, comparing them to analytical ing discussions on numerical methods and turbulence. predictions and experimental measurements. The US Department of Energy/Basic Energy Sciences Ourprimarytheoreticalfocusinourplasticitystudy[3] grantno. DE-FG02-07ER46393supportedourwork,and has been to analyze power-law correlation functions for theUSNationalCenterforSupercomputingApplications the dislocation density, plastic distortion tensor, and lo- partially provided computational resources on the Lin- cal crystalline orientation. As Fig. 5 shows, like tur- coln and Abe clusters through grant no. MSS090037. bulence simulations, these statistical properties seem to converge nicely in the continuum limit. It’s worth noting that, in both cases, non-convergence emergeswhensmallscalefeaturesappearonthewall(see ∗ Electronic address: [email protected] Fig. 6) or the interface (see Fig. 7): † Electronic address: [email protected] In the case of our simulations of plastic flow (Eq.1), [1] T. Mura, Micromechanics of defects in solids (Springer, starting from smooth initial density profiles, finite time 1987), 2nd ed. singularities develop in the form of δ-shocks. The ex- [2] P. Hähner, K. Bay, and M. Zaiser, Phys. Rev. Lett. 81, istence of finite-time singularities was shown in a 1D 2470 (1998). [3] Y.S.Chen,W.Choi,S.Papanikolaou,andJ.P.Sethna, variant of these equations, which is associated with the Phys. Rev.Lett. 105, 105501 (2010). Burgers equation [11]. Figure 6 shows how this effect [4] S. Limkumnerd and J. P. Sethna, Phys. Rev. Lett. 96, occurs by considering the two-norm differences (the in- 095503 (2006). tegrand in space of the L2 norm discussed earlier). At [5] A. Acharya,J. Mech. Phys.Solids 49, 761 (2001). t = 0.04 (see Fig. 2a) when the N = 512 curve starts to [6] A. Kurganov, S. Noelle, and G. Petrova, SIAM J. Sci. cross all the other curves, singular features start to ap- Comput. 23, 707 (2002). pear around a wall (Fig. 6b). Although the boundaries [7] A. Mignone, G. Bodo, S. Massaglia, T. Matsakos, O. Tesileanu, C. Zanni, and A. Ferrari, Astrophys. J. arenon-convergentwhenspecific locationsandtimesare Suppl.S. 170, 228 (2007). considered (Fig. 2a), the statistical properties and asso- [8] D. C. Tan, T. Zhang, T. Chang, and Y. X. Zheng, J. ciated self-similarity (Fig. 5) are convergent. Differ. Equations 112, 1 (1994). In the case of Navier-Stokes simulations (see Fig. 7), [9] J. Shi, Y.-T. Zhang, and C.-W. Shu, J. Comput. Phys. 6 (a)h=1/256 (b)h=1/512 (c)h=1/256 (d)h=1/512 t=0.04 t=1.00 FIG. 6: Non-convergence and singularity for CDD For (a) and (b), t=0.04; for (c) and (d), t=1.00. The two-norm difference between h and h/2 are plotted. At short time, t=0.04, the differences are small: (a) is empty and (b) nearly so. At later times, t=1.00, the two-norm difference becomes significant esepcially where the walls are forming (see Fig. 3 for wall locations). (a)h=1/64 (b)h=1/128 (c)h=1/256 (d)h=1/64 (e)h=1/128 (f)h=1/256 t=2.0 t=3.0 FIG. 7: Non-convergence and vortices for Navier-Stokes For (a), (b), and (c) t=2.0; for (d), (e), and (f) t=3.0, corresponding to the two red arrows in Figure 2b. The two-norm difference between h and h/2 is plotted. At t=2.0, small structures start to become strong,as in (c), and the same is true at t=3.0 for h=1/128,as in (e). 186, 690 (2003). [12] C.BardosandD.Lannes,e-printarXiv:1005.5329(2010). [10] U. Frisch, Turbulence : the legacy of A.N. Kolmogorov [13] A. Pumir and E. D. Siggia, Phys. Rev. Lett. 68, 1511 (Cambridge University Press, 1995). (1992). [11] S. Limkumnerd and J. P. Sethna, J. Mech. Phys. Solids 56, 1450 (2008).