Shock and Rarefaction Waves in Generalized Hertzian Contact Models
Hiromi Yasuda, Christopher Chong, Jinkyu Yang, Panayotis G. Kevrekidis
SShock and Rarefaction Waves in Generalized Hertzian Contact Models
H. Yasuda, C. Chong, J. Yang, and P. G. Kevrekidis Aeronautics & Astronautics, University of Washington, Seattle, WA 98195-2400, USA Department of Mathematics, Bowdoin College, Brunswick, ME 04011, USA Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003-4515, USA (Dated: November 7, 2018)In the present work motivated by generalized forms of the Hertzian dynamics associated with granular crys-tals, we consider the possibility of such models to give rise to both shock and rarefaction waves. Dependingon the value p of the nonlinearity exponent, we find that both of these possibilities are realizable. We use aquasi-continuum approximation of a generalized inviscid Burgers model in order to predict the solution profileup to times near the shock formation, as well as to estimate when it will occur. Beyond that time threshold,oscillations associated with the discrete nature of the underlying model emerge that cannot be captured bythe quasi-continuum approximation. Our analytical characterization of the above features is complemented bysystematic numerical computations. PACS numbers: 45.70.-n 05.45.-a 46.40.Cd
I. INTRODUCTION
Over the last two decades, the examination of granular crys-tals has received considerable attention, as is now summa-rized in a wide range of reviews [1–5]. Granular crystalsconsist of closely packed arrays of particles typically mod-eled as interacting with each other elastically via so-calledHertzian contacts. The resulting force depends on the geom-etry of the particles, the contact angle, and elastic propertiesof the particles [6]. Part of the reason for the wide appealof these systems is associated with their remarkable tunabil-ity, which permits to access weakly (via application of a pre-stress/precompression and for strains considerably lower thanthe precompression) or strongly (in the absence of such pre-compression) nonlinear regimes. At the same time, it is possi-ble to easily access and arrange the media in homogeneous orheterogeneous configurations. Notable examples of the laterinclude periodically arranged chains or those involving disor-der. It is, thus, not surprising that granular crystals have beenexplored as a prototypical playground for a variety of appli-cations including, among others, shock and energy absorb-ing layers [7–10], actuating devices [11], acoustic lenses [12],acoustic diodes [13] and switches [14], and sound scramblers[15, 16].Much attention has been paid to special nonlinear solutionsof the granular chain. Perhaps the most prototypical exampleexplored in the context of granular crystals (with or withoutprecompression) is the traveling solitary wave [1, 2]. Morerecently, both theoretical/computational and experimental in-terest has been expanded to another broad class of solutions,so-called bright and dark discrete breathers, which are expo-nentially localized in space and periodic in time [3, 4].Our aim in the present work is to revisit a far less exploredfamily of waveforms, namely shock waves. Shock waves arecharacterized by an abrupt, nearly discontinuous change inthe wave [17]. In the context of partial differential equations(PDEs), such as the inviscid Burgers equations, this defini-tion can be made more precise to be a solution that developsinfinite derivatives in finite time [18]. Shock-like structureshave been studied experimentally in homogeneous granular chains in the absence of precompression [19, 20]. In thoseworks, a shock wave was generated by applying a velocity toa single particle [19] or by imparting velocity continuouslyto the end of a chain [20]. More generally, in the contextof the celebrated Fermi-Pasta-Ulam (FPU) lattices (of whichthe granular chain is a specific example of), dispersive shockwaves were examined numerically for the case of general con-vex FPU potentials for arbitrary Riemann (i.e., jump) initialdata [21].In this study, we generalize the approach taken in the workof [22]. There, it was recognized that the quasi-continuumform of the Hertzian model directly relates to a second-orderPDE; for a rigorous justification, see [23]. Considering thenthe first-order nonlinear transport PDEs that represent theright and left moving waves, one retrieves effective modelsof the generalized family of the inviscid Burgers type [24].This enables the use of characteristics in order to predict theevolution of initial data, as well as the potential formation ofshock waves.The most commonly studied granular crystals are thoseconsisting of strain-hardening materials. This implies that thenonlinear exponent satisfies p > in F ∝ δ p , where F and δ are compressive force and displacement. For example, in theHertzian case of spherical particles, the nonlinear exponent is p = 3 / [25]. More recently, generalizations of the Hertziancontact law have been explored in the context of mechanicalmetamaterials, including those with strain-softening behavior(where < p < ) [26]. Examples of this include tenseg-rity structures [27] and origami metamaterial lattices [28, 29].Motivated by those recent works, we will examine the for-mation of shock waves for general values of p . We find afundamentally different dynamical behavior for the two casesof p > and p < . For p > where, for monotonicallydecreasing initial data, a shock forms due to the larger am-plitudes traveling faster from the smaller ones. In the caseof p < , the same initial data lead to a rarefaction waveform,where parts of the wave with small amplitude travel faster thanlarger amplitude parts of the wave [30]. On time scales wherethe quasi-continuum approximation remains valid, we are ableto analytically follow the solutions for both strain-hardening a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec ( p > ) and strain-softening materials ( p < ). In the vicinityof the shock time (which we predict in reasonable agreementwith the numerics from the quasi-continuum approximation),the discrete nature of the underlying lattice emerges and sig-nificantly affects the dynamics. Our systematic simulationsillustrate – by varying quantities like the precompression orthe nonlinearity power – how these predictions become pro-greessively less acurate as the model deviates from its linearanalogue.Our presentation of the above results will be structured asfollows. In section II, we will provide the theoretical back-ground and analytical findings associated with the generalizedmodel. In section III, we will compare theoretical predictionsto direct numerical simulations of the particle model. Finally,in section IV, we summarize our findings and present somepossible directions for future work. II. THEORETICAL ANALYSIS
Let u n be the relative displacement from equilibrium of the n -th particle in the lattice. As discussed in the introduction,we are motivated by the generalization of the Hertzian contactmodel to arbitrary powers not only due to applications [27–29], but also due to significant theoretical developments (e.g.,regarding traveling waves) [22, 31, 32]. In that light, we con-sider the following equations of motion: ¨ u n = [ δ + u n − − u n ] p + − [ δ + u n − u n +1 ] p + , (1)where δ is the precompression (static load) applied to thesystem and the overdot represents differentiation with respectto normalized time. Defining the strain as y n = u n − − u n ,we rewrite Eq. (1) as ¨ y n = ( δ + y n − ) p + ( δ + y n +1 ) p − δ + y n ) p . (2)Following [23] towards the derivation of a long-wavelengthapproximation characterized by smallness parameter ( ε ) andassociated spatial and temporal scales, respectively, X = εn and T = εt , we express the strain as y n ( t ) = Y ( X, T ) andobtain: ε ∂ T Y = { δ + Y ( X + ε ) } p + { δ + Y ( X − ε ) } p − { δ + Y ( X ) } p . Hence, the continuum model follows: ∂ T Y = ∂ X { ( δ + Y ) p } . (3)Now, in the spirit of [24], consider the first order (nonlineartransport) PDEs that would be compatible with Eq. (3). ∂ T Y ± α∂ X { ( δ + Y ) c } = 0 , (4)where ± indicates the two propagation directions. The param-eters α and c will be chosen such that solutions of Eq. (4) aresolutions of Eq. (3). From this we infer: ∂ T Y = ∂ T [ ∓ α∂ X { ( δ + Y ) c } ] = ∓ α c c − ∂ X (cid:110) ( δ + Y ) c − (cid:111) . (5) Comparing Eq. (3) and the right hand side of the last equalityof Eq. (5), we obtain c = ( p + 1) / and α = c − c . Usingthese definitions of c and α , we re-write Eq. (4) as ∂ T Y ± √ p ( δ + Y ) p − ∂ X Y = 0 . (6)Using then the standard technique of characteristics (see,e.g., [30] for an elementary discussion) in order to solveEq. (6), we obtain the equation along characteristic lines(along which the solution is constant) of the form: dXdT = φ ( Y ( X, T )) (7)where φ ( Y ( X, T )) = √ p ( δ + Y ( X, T )) p − . (8)Since solutions of Eq. (6) are constant along characteristiclines, Eq. (7) becomes dXdT = X − X T = φ ( Y ( X, T )) = φ ( Y ( X , . (9)In this study, we choose Y ( X ) = a sech ( bX ) (10)as the initial function, motivated by the interest in localizedinitial data.As is well known in this broad class of generalized invis-cid Burgers models, despite the smooth initial distribution ofthe strain ( y n ), a shock wave (a wave breaking effect) can beobserved to form in a finite time; i.e., the solution develops in-finite derivatives in a finite time in the continuum limit of theproblem [18]. This shock wave is created (due to the resultingmulti-valuedness) when the characteristic lines intersect. Topredict when the shock wave is formed, i.e., the shock wavetime ( T s ), we consider the following two (arbitrary) charac-teristic lines. X ( T ) = φ ( X ) T + X (11) X ( T ) = φ ( X ) T + X (12)where X , X ∈ R and X = X + h . When these two linesintersect, we obtain T = − X − X φ ( X ) − φ ( X ) = − hφ ( X + h ) − φ ( X ) . (13)The shock wave is formed at the time at which the characteris-tic lines intersect for the first time. Therefore, the shock timeis calculated as follows: T s = min (cid:20) − hφ ( X + h ) − φ ( X ) (cid:21) = 1min (cid:104) − φ ( X + h ) − φ ( X ) h (cid:105) . Then, considering h → , we obtain T s = − (cid:104) dφ ( X ) dX (cid:105) . (14)In this study, our semi-analytical prediction based on theabove considerations consists of evaluating dφ ( X ) dX and (nu-merically) identifying the relevant minimuum, which leads viaEq. (14) to a concrete estimate of T s to be compared with di-rect numerical simulations. III. NUMERICAL COMPUTATIONS AND COMPARISON
Armed with the above theoretical considerations, we nowturn to a comparison of the discrete model and the continuummodel in terms of the formation of shock waves (as a wayof quantifying the accuracy of our prediction). In order toinitialize the direct simulation of the discrete model, the initialvelocity of each particle is needed. Based on the initial strain Y ( X ) , the initial velocity can be computed as (cid:20) dy n dt (cid:21) t =0 = ddt [ Y ( X, T )] t =0 = ± ε (cid:104) √ p ( δ + Y ) p − ∂ X Y (cid:105) . The initial strain is given by Eq. (10) with the numerical con-stants a = 0 . and b = 0 . . These simulations will be com-pared to the continuum model approximations. The spatialprofile of the continuum approximation at a given time t canbe computed by solving the following implicit equation, Y ( εn, εt ) = Y ( εn − φεt ) , (15)where the wave velocity φ = ± (cid:112) p ( δ + Y ( εn, εt )) p − fol-lows directly from Eq. (6).Figure 1 shows the strain wave propagation for the standardcase of spherical Hertzian contacts p = 1 . [1, 2], while Fig. 2motivated by the works of [27–29] shows the case of p = 0 . .In these two figures, the left panels show the wave shape atfour different time instances, and the right panels show thespace-time contour plots of strain wave propagation. In eachcase the top row represents the case of no-precompression(fully nonlinear case), while progressively the middle and bot-tom row introduce higher precompression rendering the prob-lem progressively more nonlinear. For p = 1 . , the wavebreaks from its front part (see Fig. 1), whereas the p = 0 . case shows the wave breaking from the tail part as shown inFig. 2. This can be understood since the speed of propagation φ ( Y ( X, T )) in the system depends on the amplitude of thewave, see Eq. (8). For p > , points along the initial conditionwith larger amplitudes travel faster than points with smalleramplitudes. Thus, the large amplitude part of the wave over-takes the smaller amplitude part, leading to a shock formationin the front (monotonically decreasing part of the data), and ararefaction in the back (monotonically increasing part of thedata). On the contrary, for p < , the situation gets reversedas is natural to infer from the speed expression φ ( Y ( X, T )) :the rarefaction forms in the front, while the wave breakingemerges in the back.In Figs. 1 and 2, in addition to showing the dynamics ofthe original underlying lattice of Eq. (1), the prediction basedon the theoretical analysis of the nonlinear transport PDE is FIG. 1: Strain wave propagation for p = 1 . and (a) δ = 0 , (b) (c) shown, see Eq. (15). It can be inferred that the two closelymatch until the time where a discontinuity is about to format which time the discrete model starts to feature oscillatorydynamics, a dispersive feature absent in the continuum, longwavelength model. A specific diagnostic in that connectionthat we use in order to compare the discrete and continuumcases is the shock formation time T s . To obtain this time fromthe simulations, we examine the slope of the wave shape bycalculating the sign of the differences of the strains betweenadjacent particles, and we define the shock wave time ( T S )when the sign changes multiple times, i.e., when the discretecharacter of the model prevents (the continuum) shock waveformation. The resulting comparison of T s is a principal quan-titative finding of the present work.Figure 3 shows the comparison between numerical and ana-lytical shock wave time for different exponent values ( p ). Thedifference increases if p approaches unity because p = 1 in-dicates that the system becomes a linear advection equationwhich does not support shock waves. It is for that reason that FIG. 2: Strain wave propagation for p = 0 . and (a) δ = 0 , (b) (c) the shock formation time diverges as the limit p → is ap-proached. Nevertheless, the quantitative trend between thetwo cases ( p > and p < ) is well captured in both pan-els of the figure, i.e., for different values of the nonlinearity.In addition, we analyze the effect of pre-compression ( δ ) onthe shock wave time as shown in Fig. 4. As we increase thepre-compression (i.e., the system enters the weakly nonlin-ear regime and progressively approaches linear regime), thedifference between numerical and analytical results increases.One can argue that this disparity may be partially related to theway we quantify T s . In particular, as the system approaches p = 1 or δ large, linear modes become progressively moreaccessible to it (a situation to be contrasted with the sonicvacuum where no such modes exist). As a result, while in thesonic vacuum – in the effective absence of linear modes – thetime of the shock nearly coincides with the emergence of non-monotonicity near the wave front, in the progressively nearerlinear case, the non-monotonicity emerges earlier. Hence thenumerical evaluation of T s (based on the emergence of dis- creteness induced oscillations) suggests that T s is calculatedto be well before the wave breaking event predicted by the an-alytics. Despite this distinction, it is fair to say that the quasi-continuum theory provides a very good tool for predicting thepre-shock wavefront evolution and for estimating the shocktime, except in the vicinity of linear dynamical evolution. (a) δ = 0.001 (b) δ = 0.01 FIG. 3: Effect of the exponent value ( p ) on the shock wave time for (a) δ = 0 . and (b) δ = 0 . . The black square indicates theshock wave time from numerical simulations ( T S ), and the red circleis the analytical prediction ( T s ). (a) p = 0.5 (b) p = 1.5 FIG. 4: Effect of the pre-compression ( δ ) on the shock wave timefor (a) p = 0 . and (b) p = 1 . . The black square indicates theshock wave time from numerical simulations ( T S ), and the red circleis the analytical prediction ( T s ). Data is plotted on log-log scale. IV. CONCLUSIONS & FUTURE CHALLENGES
In the present work, we have investigated shock and rar-efaction wave formation in generalized granular crystal sys-tems with Hertzian contacts. Motivated by recent develop-ments in origami metamaterials [28, 29] and tensegrity struc-tures [27], we explored both the scenario of p > (including p = 3 / of spherical contacts), and that of p < of strain-softening media. We identified the emergence of both shockwaves and rarefaction waves and found a complementaritybetween the two cases. In the strain-hardening case, shockwaves arise out of monotonically decreasing initial conditionswhile rarefactions out of monotonically increasing ones. Thereverse occurs in the case of strain-softening p < me-dia. The generalization of the quasi-continuum formulationof [23] in the spirit of the nonlinear transport equations pro-posed by [24] provided us with an analytical handle in orderto characterize the evolution of pulse-like data in the strainvariables. A key quantitative prediction concerned the time offormation of the discontinuity as characterized at the quasi-continuum level. We discussed how the discreteness of thesystem, once visiting scales comparable to the lattice spacing, intervenes by generating oscillations and departing in this wayfrom the quasi-continuum description.This work suggests a number of exciting possibilities forthe future. One interesting avenue is that of obtaining infor-mation about dispersive shock waves either at the level of theKdV model, or at that of the Toda lattice and then. Then, inthe spirit of [33], we can use these to characterize the forma-tion of dispersive shock waves that emerge in the presence ofprecompression. That being said, it is unclear what becomesof such dispersive shock waves when precompression is weakor absent, as especially in the latter case there is no definitivecharacterization of dispersive shock waves. However, the di-rection of the first order equations could be an especially prof-itable one in that context. The work of [34] (see also numer-ous important references therein) suggests that for the invis-cid Burgers problem shock waves have been extensively stud-ied both analytically and numerically. It is conceivable that a (judiciously selected) discretization of the first order PDEsconsidered herein could offer considerable insight in the for-mation of shock waves, both in the weak and even in the caseof vanishing precompression limit. These are important direc-tions that are under current consideration and will be reportedin future publications. Acknowledgments
H.Y. and J.Y. acknowledge the support of the NSF underGrant No. 1553202. The work of C.C. was partially fundedby the NSF under Grant No. DMS-1615037. P.G.K grate-fully acknowledges the support of the ERC under FP7; MarieCurie Actions, People, International Research Staff ExchangeScheme (IRSES-605096), and the Alexander von HumboldtFoundation. J.Y. and P.G.K. are grateful for the support of theARO under grant W911NF-15-1-0604. [1] V. F. Nesterenko,
Dynamics of Heterogeneous Materi-als (Springer-Verlag, New York, NY, USA, 2001), ISBN1567205437.[2] S. Sen, J. Hong, J. Bang, E. Avalos, , and R. Doney, Phys. Rep. , 21 (2008).[3] G. Theocharis, N. Boechler, and C. Daraio, in
Acoustic Meta-materials and Phononic Crystals (Springer-Verlag, Berlin, Ger-many, 2013), pp. 217–251.[4] M. A. Porter, P. G. Kevrekidis, and C. Daraio, Physics Today , 44 (2015).[5] A. F. Vakakis, in Wave Propagation in Linear and NonlinearPeriodic Media (International Center for Mechanical Sciences(CISM) Courses and Lectures) (Springer-Verlag, Berlin, Ger-many, 2012), p. 257.[6] K. L. Johnson,
Contact Mechanics (Cambridge UniversityPress, Cambridge, UK, 1985).[7] C. Daraio, V. F. Nesterenko, E. B. Herbold, and S. Jin, Phys.Rev. Lett. , 058002 (2006).[8] J. Hong, Phys. Rev. Lett. , 108001 (2005).[9] F. Fraternali, M. A. Porter, and C. Daraio, Mech. Adv. Mat.Struct. , 1 (2010).[10] R. Doney and S. Sen, Phys. Rev. Lett. , 155502 (2006).[11] D. Khatri, C. Daraio, and P. Rizzo, in Nondestructive Char-acterization for Composite Materials, Aerospace Engineering,Civil Infrastructure, and Homeland Security 2008 (2008), proc.SPIE 6934, 69340U.[12] A. Spadoni and C. Daraio, Proc. Nat. Acad. Sci. USA , 7230(2010).[13] N. Boechler, G. Theocharis, and C. Daraio, Nature Materials , 665 (2010).[14] F. Li, P. Anzel, J. Yang, P. G. Kevrekidis, and C. Daraio, Nat.Comm. , 5311 (2014).[15] C. Daraio, V. F. Nesterenko, and S. Jin, Phys. Rev. E , 016603(2005).[16] V. F. Nesterenko, C. Daraio, E. B. Herbold, and S. Jin, Phys.Rev. Lett. , 158702 (2005).[17] M. J. Ablowitz and M. Hoefer, Scholarpedia , 5562 (2009).[18] J. Smoller, Shock Waves and Reaction–Diffusion Equations (Springer-Verlag, Berlin, Germany, 1983).[19] E. B. Herbold and V. F. Nesterenko, Phys. Rev. E , 021304(2007).[20] A. Molinari and C. Daraio, Phys. Rev. E , 056602 (2009).[21] M. Herrmann and J. D. M. Rademacher, Nonlinearity , 277(2010).[22] K. Ahnert and A. Pikovsky, Phys. Rev. E , 026209 (2009).[23] C. Chong, P. G. Kevrekidis, and G. Schneider, Disc. Cont. Dyn.Sys. A , 3403 (2014).[24] B. E. McDonald and D. Calvo, Phys. Rev. E , 066602 (2012).[25] H. Hertz, J. Reine. Angew. Math , 156 (1881).[26] E. B. Herbold and V. F. Nesterenko, Phys. Rev. Lett. ,144101 (2013), URL http://link.aps.org/doi/10.1103/PhysRevLett.110.144101 .[27] F. Fraternali, G. Carpentieri, A. Amendola, R. E. Skelton,and V. F. Nesterenko, Applied Physics Letters , 201903(2014), URL http://scitation.aip.org/content/aip/journal/apl/105/20/10.1063/1.4902071 .[28] H. Yasuda and J. Yang, Phys. Rev. Lett. , 185502(2015), URL http://link.aps.org/doi/10.1103/PhysRevLett.114.185502 .[29] H. Yasuda, C. Chong, E. G. Charalampidis, P. G. Kevrekidis,and J. Yang, Phys. Rev. E , 043004 (2016).[30] W. Strauss, Partial Differential Equations: An Introduction (John Wiley & Sons, Hoboken, 2008).[31] G. James and D. Pelinovsky, Proceedings of the RoyalSociety of London A: Mathematical, Physical andEngineering Sciences (2014), ISSN 1364-5021,http://rspa.royalsocietypublishing.org/content/470/2165/20130462.full.pdf,URL http://rspa.royalsocietypublishing.org/content/470/2165/20130462 .[32] E. Dumas and D. E. Pelinovsky, SIAM J. Math. Anal. , 4075(2014).[33] Y. Shen, P. G. Kevrekidis, S. Sen, and A. Hoffman, Phys. Rev.E , 022905 (2014).[34] C. V. Turner and R. R. Rosales, Stud. Appl. Math.99