Soft materials with engineered microstructure support nonlinear waves which can be har- nessed for various applications, from signal communication to impact mitigation. Such waves are governed by nonlinear coupled diﬀerential equations whose analytical solution is seldom trackable, hence emerges the need for suitable numerical solvers. Based on a ﬁnite-volume method in one space dimension, we here develop a designated scheme for nonlinear waves with two coupled components that propagate in soft laminates. We apply our scheme to a periodic laminate made of two alternating compressible Gent layers, and consider two cases. In one case, we analyze a motion whose component along the lamination direction is coupled to a component in the layers plane, and discover vector solitary waves in a continuum medium. In the second case, we analyze a motion with two coupled components in the plane of the layers, and observe a train of linearly polarized solitary waves, followed by a single circularly polarized wave. The framework we developed oﬀers a platform for further investigation of these waves and their extension to higher dimensional problems. Keywords : ﬁnite-volume method, composites, laminates, soft materials, nonlinear elasticity, ﬁnite deformations, wave propagation, solitary waves, shocks, Gent material 1 Introduction Highly deformable materials with architectured microstructure are nowadays accessible by virtue of the current manufacturing abilities ( Bandyopadhyay et al. ,2015,Truby and Lewis ,2016). These materials exhibit geometrical and constitutive nonlinearities that give rise to rich physics ( Bertoldi et al.,2017), and in particular diverse wave phenomena such as unidirectional transmission, shocks ∗Corresponding author. Tel.: +1 972 778871613. E-mail address : meshmuel@technion.ac.il (G. Shmuel). 1 and self-reinforcing waves ( Nadkarni et al. ,2014,2016,Samsonov et al. ,2017,Ziv and Shmuel , 2019,Deng et al. ,2019a ,Katz and Givli ,2019). Waves in nonlinear solids have recently regained scientiﬁc and technological interest, owing to the realization that their features can be harnessed for various applications, such as signal transmission ( Raney et al. ,2016), impact mitigation ( Yasuda et al.,2019), energy harvesting ( Hwang and Arrieta ,2018) and mechanical diodes ( Deng et al. , 2018). The mathematical modeling of waves in nonlinear solids with microstructure is given by nonlin- ear coupled partial diﬀerential equations, whose analytical solution is seldom trackable ( Davison , 1966,Carroll ,1967,Destrade and Saccomandi ,2005,Chockalingam and Cohen ,2020), hence the need for designated numerical solvers. Notably, these solvers should account for the fact that when the deformations are large there is a distinction between the reference and current conﬁguration space, and therefore employ the so-called Lagrangian formulation ( Drumheller ,1998,Kluth and De- sprés,2010,Bonet et al. ,2015). Among the diﬀerent schemes developed for nonlinear Lagrangian elastodynamics, we recall those based on the ﬁnite diﬀerence method ( Aboudi and Benveniste , 1973), ﬁnite element method ( Bou Matar et al. ,2012), and speciﬁcally the ﬁnite-volume method (Berjamin et al. ,2018). Finite-volume methods are useful in solving problems whose physics is governed by space- time conservation laws ( Godlewski and Raviart ,1996,Toro,1997,Trangenstein ,2009), and as such are useful in elastodynamics ( LeVeque ,2002a ). Importantly, these methods are conservative in the sense that by construction they satisfy the conservation law in an exact manner—a useful feature in problems whose solution is discontinuous ( Lax and Wendroﬀ ,1960). Based on the solution of a Riemann problem at the interface between grid cells, LeVeque (1997) developed a general high-resolution scheme for nonlinear hyperbolic systems which depend on multiple space dimensions and are not in conservation form. Such a deviation from the conservation form can occur when the medium is heterogeneous—the type of medium we focus on in this work, giving rise to a spatially-varying ﬂux. The application of ﬁnite-volume algorithms based on solutions of Riemann problems has been extended to thermoelastic materials ( Berezovski and Maugin ,2001), elastoplastic materials ( Barton et al. ,2010), and materials that exhibit softening ( Berjamin et al. , 2018), among other media. Fogarty and LeVeque (1999) have demonstrated the eﬃciency of Leveque’s approach in diﬀer- ent problems of acoustic waves and one- and two-dimensional heterogeneous linear media. The term acoustic denotes waves in media that cannot sustain shear, hence corresponds to pressure (compression) waves alone. In the periodic case, they have also reﬁned the method using a new limiter function —a scalar function aimed at limiting the gradient of the solution near discontinu- ities—relatively to the one used by LeVeque (1997). Later on, LeVeque (2002b ) and Bale et al. (2002) have adapted the approach to pressure waves in one- and two-dimensional nonlinear hetero- 2 geneous media, by decomposing the ﬂux diﬀerence at the interface between grid cells. This is in contrast with the standard approach that decomposes the conserved vector. These works revealed solitary waves, namely, nonlinear waves that are able to maintain constant speed and proﬁle by virtue of a balance between nonlinear and dispersive features of the system ( Dauxois and Peyrard , 2006,Hussein and Khajehtourian ,2018). These one-dimensional elastic solitary waves were later studied in greater detail by LeVeque and Yong (2003a ,b); see also the works of Engelbrecht et al. (2007),Xu et al. (2007). Here, we extend the ﬂux-based wave decomposition scheme of LeVeque (2002b ) and Bale et al. (2002) to nonlinear elastic heterogeneous media exhibiting coupled ﬁnite-amplitude waves. In contrast with acoustics where there is only volumetric waves and the only stress component is in the propagation direction, in elastodynamics there are also shear waves which may couple with volumetric waves, and yield traction components in the plane transverse to the propagation direction ( Lustig et al. ,2019). The media of interest are composites made of soft layers, where the component of the deformation in the lamination direction (referred to as the axial component) is coupled with a transverse component in the plane of the layers ( Davison ,1966). We also consider a second case in which two components in the layers plane are coupled one with the other ( Collins , 1966),i.e., the coupling of two shear waves. By developing of a suitable matrix representation, our extension accounts for the generation of such additional waves and their coupling. To the best of our knowledge, while ﬁnite-volume schemes have been developed for two-dimensional nonlinear elastodynamics, e.g., byBerjamin et al. (2019), the case of spatially-varying ﬂux in that setting has yet to be explicitly addressed 1. We validate our method using two test cases. In the ﬁrst test, we consider one layer bounded between two semi-inﬁnite layers with diﬀerent parameters, and note that these layers respond nonlinearly to ﬁnite-amplitude deformations. We subject the middle layer to an initial small- amplitude shear strain and show that the numerical solution captures the response as predicted by the analytical solution of the limiting linear problem ( Shmuel and Moiseyev ,2020). The second test is of nonlinear waves and hence more challenging, where we address the canonical problem of wave scattering at an interface between two nonlinearly elastic half-spaces; here, however, the waves are of ﬁnite-amplitude and the displacements are coupled. The aspect of shock evaluation in this setting was addressed by Davison (2008), who studied axial waves. Here, we extend the settings to shocks associated with waves comprising two coupled components, noting that to the best of our knowledge this is the ﬁrst numerical experiment of such problems. An analytical solution to this problem is not available, and hence our ﬁnite-volume based solution to the linearized algebraic equations is compared with numerical solutions using Newton’s method to the exact nonlinear 1While the setting of Berjamin et al. (2019) is two-dimensional nonlinear elastodynamics, the numerical method they provide when the coeﬃcients are spatially-varying is for the one-dimensional case (Appendix B therein). 3 equations, to ﬁnd an excellent agreement. The model we use to describe the nonlinear response of the media in our test cases is the compressible Gent model ( Gent,1996). This model—originally developed for capturing the strain stiﬀening of rubber at large strains—was recently shown useful by Ziv and Shmuel (2019) for modeling shear shocks that were experimentally observed by Catheline et al. (2003),Espíndola et al. (2017), and modeling tensile-induced shocks that were experimentally observed by Niemczura and Ravi-Chandar (2011). Subsequently, we apply our scheme to a pre-strained laminate made of two alternating com- pressible Gent layers. Remarkably, in the case where the axial and transverse components of the motion are coupled, our numerical solution reveals a generation of elastic vector solitary waves. As mentioned, solitary waves are nonlinear waves that propagate with constant speed and proﬁle by virtue of a balance between nonlinearity and dispersion in the system. The term vector refers to the case where these solitary waves consist of two (or more) components andpolarizations that are coupled one with the other ( Chen ,1994,Chen et al. ,2012). To the best of our knowledge, our results are the ﬁrst report of vector solitary waves in an elastic continuum , based on the equations of nonlinear elastodynamics. Our observation is preceded by the ﬁrst construction of vector solitary waves in discrete mechanical systems that were conceived by Deng et al. (2017,2019b ). There, the model is a periodic repetition of rigid squares that are interlinked by springs, thereby supporting transitional and rotational waves. Interestingly, while these mechanical waves are slower at higher amplitudes, the vector solitary waves in our continuum laminated medium are faster at higher amplitudes, similarly to the KdV solitons and the solitary waves analyzed by LeVeque and Yong (2003b ), as we show in what follows. In the second case where the coupling is between two transverse components, we discover the propagation of (scalar) solitary waves that are linearly polarized, namely, their polarization direction is ﬁxed. We also observe the propagation of a single slower wave which is circularly polarized, namely, its direction rotates in the transverse plane. Our results are presented in the following order. Sec. 2summarizes the equations governing motions with two coupled components in laminates made of compressible Gent layers, together with the derivation of the Gent phase velocities in each case. Sec. 3details our ﬁnite-volume method for the solution of these equations. The validation of our method using the two test cases is provided in Sec.4, and the numerical experiments of coupled motions in periodic Gent layers is carried out in Sec.5. We summarize our main results and conclusions in Sec. 6, together with an outlook towards future work. 4 2 Governing equations We use the standard governing equations in Lagrangian continuum mechanics to formulate the problem of interest, see, e.g., the books of Ogden (1997),Bonet and Wood (1997). We consider a soft composite made of two perfectly bonded alternating hyperelastic phases that are laminated in theX1direction. The layers are inﬁnite in the X2andX3directions; the loads in these directions are assumed to be uniform, and hence the response of the laminate is only a function of X1. We denote the initial thickness, mass density and strain energy density function of layer nbyH(n),ρ(n) L andΨ(n). In the numerical simulations to follow, we will use the compressible 2Gent energy function (Gent,1996,Lopez-Pamies and Castañeda ,2007) to model the phases, given by Ψ(n)(F)=−µ(n)J(n) m 2ln/parenlefttpA /parenleftbtA1−trFTF−3 J(n) m/parenrighttpA /parenrightbtA−µ(n)lndet F+/parenlefttpA /parenleftbtAκ(n) 2−µ(n) 3−µ(n) J(n) m/parenrighttpA /parenrightbtA(detF−1)2,(1) where Fis the deformation gradient (to be deﬁned formally later), µ(n)andκ(n)correspond to the shear and bulk moduli in the limit of small strains, respectively, and J(n) mmodels the strain stiﬀening that rubber exhibits due to the limited extensibility of its polymer chains. In the limit Jm→∞, the Gent model reduces to the fundamental neo-Hookean model ( Horgan ,2015). The motion of the composite is described in terms of χ: a continuous and twice diﬀerentiable function (except at material interfaces), which maps material points from their initial position Xto their current position xat time t, such that x=χ(X,t). As mentioned, we focus on deformations that are only functions of X1, namely, x1=X1+u1(X1,t),x2=X2+u2(X1,t),x3=X3+u3(X1,t), (2) where uiare the components of the displacement ﬁeld. The problem thus amounts to determining ui(X1,t)for a given set of boundary and initial conditions. This requires the satisfaction of balance laws for the corresponding ﬁrst Piola-Kirchhoﬀ stress tensor P=∇FΨ, where F=∇Xχis the deformation gradient. For the compressible Gent model, the F-gradient yields P=µ/parenleftBigg 1−trFTF−3 Jm/parenrightBigg−1 F+/parenleftBigg κ−2µ Jm−2µ 3/parenrightBigg detF(detF−1)F−T−µF−T. (3) It is worth highlighting now what are the nonlinearities in the Cauchy stress σ=(detF)−1PFT that emerge when employing the Gent model. First, it has a term proportional to FFT, hence 2For completeness, we note that in the incompressible limit (detF≡1), the kinematic constraint is handled using the addition of a Lagrange multiplier ( Ogden ,1997,Bonet and Wood ,1997). 5 (a) (b) (c) Figure 1: (a) Reference conﬁguration of a representative phase. (b) Illustrative deformation when u1≡0, and (c)u3≡0. does not neglect quadratic terms in the components of the displacement gradient, in contrast with linear elasticity. Secondly, the coeﬃcient 3that multiplies FFTis also a nonlinear function of the deformation, in contrast with the Hookean and neo-Hookean models, where this coeﬃcient is constant. These nonlinearities play a crucial role in evolution of nonlinear waves; smooth nonlinear waves are solutions to the Lagrangian equations of motion ∇X·P=ρLχ,tt, (4) as discussed later. For the mapping of interest ( 2), Eq. ( 4) reduces to Pi1,1=ρL∂2ui ∂t2,i=1,2,3. (5) Focusing on motions in which either u1≡0 oru3≡0, we rewrite Eq. ( 5) as u1≡0 :P21,1=ρL∂2u2 ∂t2,P31,1=ρL∂2u3 ∂t2; (6a,b) u3≡0 :P11,1=ρL∂2u1 ∂t2,P21,1=ρL∂2u2 ∂t2. (6c,d) Physically, u1≡0 corresponds to a transverse motion in the plane perpendicular to X1(Fig. 1b), whereas u3≡0 corresponds to coupled axial and transverse motions (Fig. 1c). We can analyze the two cases together using a unifying equation by replacing the ﬁrst index of Pin Eqs. ( 6a) and ( 6c) byaand in Eqs. ( 6b) and ( 6d) by b, and writing 3Note that the coeﬃcient that multiplies the hydrostatic part is also a nonlinear function of the deformation. 6 Pa1,1=ρL∂2ua ∂t2,Pb1,1=ρL∂2ub ∂t2, (7) where we recall that ρLis a piecewise-constant function that varies between phases, while uiand Pi1are continuous functions 4. We denote the spatial and temporal derivatives of uiby ǫi:=∂ui ∂X1,vi:=∂ui ∂t,i=a,b, (8) and put Eq. ( 7) in the form /parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtAǫa ǫb ρLva ρLvb/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,t+/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA−va −vb −Pa1 −Pb1/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,X1=/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA0 0 0 0/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA, (9) By applying the chain rule, we obtain /parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtAǫa ǫb ρLva ρLvb/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,t+/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA0 0−1 ρL0 0 0 0−1 ρL −α−β0 0 −γ−δ0 0/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtAǫa ǫb ρLva ρLvb/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,X1=/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA0 0 0 0/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA, (10) where α=∂Pa1 ∂ǫa,β=∂Pa1 ∂ǫb,γ=∂Pb1 ∂ǫa,δ=∂Pb1 ∂ǫb, (11) and for later use, we denote the matrix in Eq. ( 10) byG. When specialized to the compressible Gent model and the case u1≡0, the Piola components required for calculating Eq. ( 11) are P21=µǫ2/parenlefttpA /parenleftbtA1−ǫ2 2+ǫ2 3 Jm/parenrighttpA /parenrightbtA−1 , P31=µǫ3/parenlefttpA /parenleftbtA1−ǫ2 2+ǫ2 3 Jm/parenrighttpA /parenrightbtA−1 ,(12) 4The continuity of Pi1results from the standard application of the balance of linear momentum in its integral form at the reference conﬁguration near material interfaces, and uiare continuous there since the layer are perfectly bonded. 7 where when u3≡0 the components are P11=µ(ǫ1+1)/parenlefttpA /parenleftbtA1−ǫ2 1+ǫ2 2+2ǫ1 Jm/parenrighttpA /parenrightbtA−1 +/parenleftBigg κ−2µ Jm−2µ 3/parenrightBigg ǫ1−µ ǫ1+1, P21=µǫ3/parenlefttpA /parenleftbtA1−ǫ2 1+ǫ2 2+2ǫ1 Jm/parenrighttpA /parenrightbtA−1 .(13) The eigenvalues of Eq. ( 10) are the following characteristic wave velocities in the material c=±/radicalBigg 1 2ρL/bracketleftBigg α+δ±/radicalBig (α−δ)2+4βγ/bracketrightBigg ≕±c±, (14) where c+andc−correspond to the plus and minus sign of the inner square root, respectively. Eq. (10) is hyperbolic when c±are real, i.e., when the conditions α+δ−/radicalBig (α−δ)2+4βγ> 0,(α−δ)2+4βγ> 0, (15) are satisﬁed ( Davison ,1966). The compressible Gent model satisﬁes conditions ( 15) for all its admissible strains before the “lock-up”/parenleftBig trFTF−3<Jm/parenrightBig , hence Eq. ( 10) is suitable for numerical solutions based on ﬁnite volume methods ( LeVeque ,2002a ). In the case of coupled axial and transverse displacements, the explicit expressions for the velocities of the compressible Gent model are c2 ±=d3 2ρL+µ 2ρLd2 1+/parenleftBig d2 1+ǫ2 2+Jm/parenrightBig µJm 2ρLd2 2±/radicalBig 48d6 1µ2ǫ2 2J2m+/parenleftBig d2 1d2 2d3+2d4 1µJm−2d2 1ǫ2 2µJm+d2 2µ/parenrightBig 2 2√ 3ρLd2 2d2 1, (16) where d1=ǫ1+1,d2=ǫ2 1+2ǫ1+ǫ2 2−Jm,d3=κ−2µ 3−2µ Jm. (17) In this case, the slower velocity c−is a monotonically increasing function of the strains, and the corresponding wave is referred to as genuinely nonlinear. The faster velocity c+has local minima at certain strains ( Ziv and Shmuel ,2019), and hence the corresponding wave is not genuinely nonlinear. In the following sections we focus on strains far from these local minima, such that this mode can be considered genuinely nonlinear. The fact that the velocities are nonlinear functions of the strain components allows the formation of shocks under certain conditions; this formation is the focus of previous works ( Ziv and Shmuel ,2019,Chockalingam and Cohen ,2020) and the relevant conditions are brieﬂy discussed in Appendix A . In the limit of linear elasticity, c−andc+associated with Eq. ( 6c,d) correspond to the velocities of shear and pressure waves, respectively. 8 Collins (1966) showed that coupled transverse deformations with ﬁnite amplitude associated with Eq. ( 6a,b) propagate in isotropic media as a combination of circularly polarized waves with the velocity c−and linearly polarized waves with the velocity c+. This nature of motion is revealed using the transformation ǫ2=ǫTcosθ,ǫ 3=ǫTsinθ, (18) which decouples the waves such that ǫTandθare constants across waves propagating with the velocities c−andc+, respectively. In other words, the circularly polarized wave is characterized by a ﬁxed strain and rotating polarization, where for the linearly polarized wave the polarization is ﬁxed and the strain magnitude varies. The corresponding velocities in a compressible Gent material are c−=/radicalBigg µJm ρL/parenleftBig Jm−ǫ2 T/parenrightBig,c+=/radicaltp/radicalvertex/radicalvertex/radicalvertex/radicalbtµJm/parenleftBig Jm+ǫ2 T/parenrightBig ρL/parenleftBig Jm−ǫ2 T/parenrightBig2, (19) respectively. Note that c−is constant sinceǫTis a constant too. Hence, circularly polarized waves are eﬀectively propagating as linear waves, and are referred to as linearly degenerate ( LeVeque , 2002a ). By contrast, the quantity ǫTvaries for linearly polarized waves, hence their velocity c+is a nonlinear function of the strain measure. Speciﬁcally, it is a monotonically increasing function of ǫT, and therefore the linearly polarized mode is genuinely nonlinear. In the subsequent sections we formulate a ﬁnite-volume method for layered materials governed by Eq. ( 10), validate it using two benchmark problems, and apply it in a numerical experiment of coupled waves in nonlinear laminates. 3 Finite-volume method First, we deﬁne the conserved vector qand its ﬂux vector faccording to q≔/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtAǫa ǫb ρLva ρLvb/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,f≔−/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtAva vb Pa1 Pb1/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA, (20) such that the two forms of the governing equations given in Eqs. ( 9)-(10) are written respectively as q,t+f,X1=0, (21a) q,t+∇qf·q,X1=0, (21b) 9 (a) f(i+1)f(i)f(i−1)q(i+1)q(i)q(i−1) ∆X ∆X ∆X (b) X1t c=−c(i) + r=r(i) 1c=−c(i) − r=r(i) 2 c=c(i+1) + r=r(i+1) 4c=c(i+1) − r=r(i+1) 3 q=q(i),f=f(i)f=f(i,i+1) q=q(i+1),f=f(i+1) Figure 2: (a) Illustration of the discretization of the material in the numerical scheme. (b) Representative characteristic curves associated with the solution to the Riemann problem of two adjacent phases between cells iandi+1. where∇qf≡G.Eq. (21b) is in the form of a conservation law for q. We consider a grid of uniform length ∆X, and approximate the values of qandfas constants within each length element using the value at the center of the element; the ithelement is denoted using superscript i(Fig. 2a). Following LeVeque (1997), we formulate a ﬁnite-volume scheme based on an approximate solution to the Riemann problem between two adjacent cells iandi+1. We adopt the ﬂux-decomposition approach of LeVeque (2002b ) and Bale et al. (2002),i.e.,we solve the system f(i+1)−f(i)=2/summationdisplay k=1/dblbracketleftAk/dblbracketrightir(i) k+4/summationdisplay k=3/dblbracketleftAk/dblbracketrightir(i+1) k; (22) here,{rk}are eigenvectors of G, namely, 10 r1=/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA1 η ρLc+ ηρLc+/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,r2=/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA−ζ 1 −ζρLc− ρLc−/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,r3=/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA−ζ 1 ζρLc− −ρLc−/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA,r4=/parenlefttpA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA/parenleftexA /parenleftbtA1 η −ρLc+ −ηρLc+/parenrighttpA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA/parenrightexA /parenrightbtA, (23) whereη=−/parenleftbigα−ρLc+/parenrightbig/β,ζ=β//parenleftbigα−ρLc−/parenrightbig, and/dblbracketleftAk/dblbracketrightidenotes the jump in fin the direction of itsktheigenvector between cells iandi+1. Thus, Eq. ( 22) constitutes a linear algebraic system of equations for/dblbracketleftAk/dblbracketrighti. Physically, this solution corresponds to four waves propagating with the velocities cn=−c+,−c−,c−andc+, which contains jumps in fandqthat are proportional to the respective eigenvectors. Characteristic curves as representative solutions to the Riemann problem of two adjacent phases in cells iandi+1 are shown in Fig. 2b, as given by Eq. ( 22). These solutions are composed of four linear waves with a discontinuity in the state and ﬂux ﬁelds, where the slopesof the characteristic lines are their velocity. These waves are separated into two rightward- and two leftward propagating waves. As noted in the introduction, the common approach is to decompose the diﬀerence q(i+1)−q(i) instead of f(i+1)−f(i). In our problem the mismatch in the material parameters between the adjacent cells results in qcontaining an additional jump at the interface which is not in the direction of any one of the eigenvectors. As advocated by LeVeque (2002b ) and Bale et al. (2002) in the context of pressure (compression) waves in nonlinear heterogeneous solids, the decomposition of the diﬀerence in fdelivers a simpler set of equations, since fmust be continuous across the interface. Next, we integrate the approximate solution in a ﬁnite-volume scheme, as done by LeVeque (1997). We use a Godunov-type method with a second order correction ( Godunov ,1959) q(i)(t+∆t)=q(i)(t)−∆t ∆X/parenleftBig f(i,i+1)−f(i−1,i)/parenrightBig −∆t 2∆X/parenleftBigˆf(i,i+1)−ˆf(i−1,i)/parenrightBig , (24) where f(i,i+1)≔f(i)+/dblbracketleftA1/dblbracketrightir(i) 1+/dblbracketleftA2/dblbracketrightir(i) 2≡f(i+1)−/dblbracketleftA4/dblbracketrightir(i+1) 4−/dblbracketleftA3/dblbracketrightir(i+1) 3(25) is the value of fat the interface (see Fig. 2b), and ˆf(i,i+1)=4/summationdisplay k=1signc(i+j) k/parenleftBigg 1−∆t ∆X/barex/barex/barex/barexc(i+j) k/barex/barex/barex/barex/parenrightBigg /dblbracketleftAk/dblbracketrightir(i+j) k(26) is the second order correction, with jequals 0 (resp. 1) for k=1,2 (resp. k=3,4). While the second order correction improves the numerical solution when the exact solution is smooth, in case of discontinuities it may lead to unwanted numerical oscillations. To resolve this 11 Setρ/bracketleftBig kg/m3/bracketrightBig µ[kPa]κ[MPa] Jm 1 1000 100 0.5 10 2 1000 400 2 10 3 500 200 1 10 4 4000 200 1 10 Table 1: Sets of Gent moduli used in the numerical simulations, chosen from the characteristic range of the values that ﬁts elastomers. issue, we employ the approach that was introduced by LeVeque (1997) which uses a wave limiter , namely, we replace /dblbracketleftAk/dblbracketrightiin Eq. ( 26) byφ/parenleftBig θ(i,i+1) k/parenrightBig /dblbracketleftAk/dblbracketrighti, whereφis a function of θ(i,i+1) k=c(i+j) k/dblbracketleftAk/dblbracketrighti−2j+1 c(i−j+1) k/dblbracketleftAk/dblbracketrightirT(i+j) k·r(i−j+1) k rT(i+j) k·r(i+j) k. (27) (See a discussion on earlier versions of wave limiters by LeVeque ,1997, and the