Observation of vector solitary waves in soft laminates using a ﬁnite-volume method

Soft materials with engineered microstructure support nonlinear waves which can be harnessed 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.


Introduction
Highly deformable materials with architectured microstructure are nowadays accessible by virtue of the current manufacturing abilities (Bandyopadhyay et al., 2015, Truby andLewis, 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 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 scientific 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 nonlinear coupled partial differential 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 configuration space, and therefore employ the so-called Lagrangian formulation (Drumheller, 1998, Kluth and Després, 2010, Bonet et al., 2015).Among the different schemes developed for nonlinear Lagrangian elastodynamics, we recall those based on the finite difference method (Aboudi and Benveniste, 1973), finite element method (Bou Matar et al., 2012), and specifically the finite-volume method (Berjamin et al., 2018).
Finite-volume methods are useful in solving problems whose physics is governed by spacetime 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 Wendroff, 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 flux.The application of finite-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 efficiency of Leveque's approach in different 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 refined the method using a new limiter function-a scalar function aimed at limiting the gradient of the solution near discontinuities-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-geneous media, by decomposing the flux difference 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 profile by virtue of a balance between nonlinear and dispersive features of the system (Dauxois andPeyrard, 2006, Hussein andKhajehtourian, 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 flux-based wave decomposition scheme of LeVeque (2002b) and Bale et al. (2002) to nonlinear elastic heterogeneous media exhibiting coupled finite-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 finite-volume schemes have been developed for two-dimensional nonlinear elastodynamics, e.g., by Berjamin et al. (2019), the case of spatially-varying flux in that setting has yet to be explicitly addressed1.
We validate our method using two test cases.In the first test, we consider one layer bounded between two semi-infinite layers with different parameters, and note that these layers respond nonlinearly to finite-amplitude deformations.We subject the middle layer to an initial smallamplitude 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 finite-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 first numerical experiment of such problems.An analytical solution to this problem is not available, and hence our finite-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 coefficients are spatially-varying is for the one-dimensional case (Appendix B therein).equations, to find 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 stiffening 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 compressible 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 profile 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 and polarizations that are coupled one with the other (Chen, 1994, Chen et al., 2012).To the best of our knowledge, our results are the first report of vector solitary waves in an elastic continuum, based on the equations of nonlinear elastodynamics.Our observation is preceded by the first construction of vector solitary waves in discrete mechanical systems that were conceived by Deng et al. (2017Deng et al. ( , 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 fixed.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. 2 summarizes 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. 3 details our finite-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.

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 the X 1 direction.The layers are infinite in the X 2 and X 3 directions; the loads in these directions are assumed to be uniform, and hence the response of the laminate is only a function of X 1 .We denote the initial thickness, mass density and strain energy density function of layer n by In the numerical simulations to follow, we will use the compressible2 Gent energy function (Gent, 1996, Lopez-Pamies andCastañeda, 2007) to model the phases, given by where F is the deformation gradient (to be defined formally later), µ (n) and κ (n) correspond to the shear and bulk moduli in the limit of small strains, respectively, and J (n) m models the strain stiffening that rubber exhibits due to the limited extensibility of its polymer chains.In the limit J m → ∞, 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 differentiable function (except at material interfaces), which maps material points from their initial position X to their current position x at time t, such that x = χ (X, t).As mentioned, we focus on deformations that are only functions of X 1 , namely, where u i are the components of the displacement field.The problem thus amounts to determining u i (X 1 , t) for a given set of boundary and initial conditions.This requires the satisfaction of balance laws for the corresponding first Piola-Kirchhoff stress tensor P = ∇ F Ψ, where F=∇ X χ is the deformation gradient.For the compressible Gent model, the F-gradient yields It is worth highlighting now what are the nonlinearities in the Cauchy stress σ = (detF) −1 PF T that emerge when employing the Gent model.First, it has a term proportional to FF T , 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 andWood, 1997).does not neglect quadratic terms in the components of the displacement gradient, in contrast with linear elasticity.Secondly, the coefficient3 that multiplies FF T is also a nonlinear function of the deformation, in contrast with the Hookean and neo-Hookean models, where this coefficient is constant.These nonlinearities play a crucial role in evolution of nonlinear waves; smooth nonlinear waves are solutions to the Lagrangian equations of motion as discussed later.For the mapping of interest (2), Eq. ( 4) reduces to Focusing on motions in which either u 1 ≡ 0 or u 3 ≡ 0, we rewrite Eq. ( 5) as Physically, u 1 ≡ 0 corresponds to a transverse motion in the plane perpendicular to X 1 (Fig. 1b), whereas u 3 ≡ 0 corresponds to coupled axial and transverse motions (Fig. 1c).We can analyze the two cases together using a unifying equation by replacing the first index of P in Eqs. ( 6a) and (6c) by a and in Eqs. ( 6b) and (6d) by b, and writing where we recall that ρ L is a piecewise-constant function that varies between phases, while u i and P i1 are continuous functions4.We denote the spatial and temporal derivatives of u i by and put Eq. ( 7) in the form By applying the chain rule, we obtain where and for later use, we denote the matrix in Eq. ( 10) by G.When specialized to the compressible Gent model and the case u 1 ≡ 0, the Piola components required for calculating Eq. ( 11) are 4The continuity of P i1 results from the standard application of the balance of linear momentum in its integral form at the reference configuration near material interfaces, and u i are continuous there since the layer are perfectly bonded.
where when u 3 ≡ 0 the components are The eigenvalues of Eq. ( 10) are the following characteristic wave velocities in the material where c + and c − 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 are satisfied (Davison, 1966).The compressible Gent model satisfies conditions (15) for all its admissible strains before the "lock-up" trF T F − 3 < J m , hence Eq. ( 10) is suitable for numerical solutions based on finite 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 where 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 andShmuel, 2019, Chockalingam andCohen, 2020) and the relevant conditions are briefly discussed in Appendix A. In the limit of linear elasticity, c − and c + associated with Eq. (6c,d) correspond to the velocities of shear and pressure waves, respectively.Collins (1966) showed that coupled transverse deformations with finite 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 which decouples the waves such that ǫ T and θ are constants across waves propagating with the velocities c − and c + , respectively.In other words, the circularly polarized wave is characterized by a fixed strain and rotating polarization, where for the linearly polarized wave the polarization is fixed and the strain magnitude varies.The corresponding velocities in a compressible Gent material are respectively.Note that c − is constant since ǫ T is a constant too.Hence, circularly polarized waves are effectively propagating as linear waves, and are referred to as linearly degenerate (LeVeque, 2002a).By contrast, the quantity ǫ T varies for linearly polarized waves, hence their velocity c + is a nonlinear function of the strain measure.Specifically, it is a monotonically increasing function of ǫ T , and therefore the linearly polarized mode is genuinely nonlinear.
In the subsequent sections we formulate a finite-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.

Finite-volume method
First, we define the conserved vector q and its flux vector f according to such that the two forms of the governing equations given in Eqs. ( 9)-( 10) are written respectively as where ∇ q f ≡ 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 q and f as constants within each length element using the value at the center of the element; the i th element is denoted using superscript i (Fig. 2a).Following LeVeque (1997), we formulate a finite-volume scheme based on an approximate solution to the Riemann problem between two adjacent cells i and i + 1.We adopt the flux-decomposition approach of LeVeque (2002b) and Bale et al. (2002), i.e., we solve the system here, {r k } are eigenvectors of G, namely, where and A k i denotes the jump in f in the direction of its k th eigenvector between cells i and i + 1.Thus, Eq. ( 22) constitutes a linear algebraic system of equations for A k i .Physically, this solution corresponds to four waves propagating with the velocities c n = −c + , −c − , c − and c + , which contains jumps in f and q that are proportional to the respective eigenvectors.Characteristic curves as representative solutions to the Riemann problem of two adjacent phases in cells i and i + 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 flux fields, 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 difference 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 q containing 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 difference in f delivers a simpler set of equations, since f must be continuous across the interface.
Next, we integrate the approximate solution in a finite-volume scheme, as done by LeVeque (1997).We use a Godunov-type method with a second order correction (Godunov, 1959) where is the value of f at the interface (see Fig. 2b), and is the second order correction, with j equals 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.issue, we employ the approach that was introduced by LeVeque (1997) which uses a wave limiter, namely, we replace A k i in Eq. ( 26) by φ θ (i,i+1) k A k i , where φ is a function of (See a discussion on earlier versions of wave limiters by LeVeque, 1997, and the references therein.) The monotonized centered function is a simple example of φ, defined by Lastly, in order for the method to be numerically stable, it is necessary for the Courant-Friedrichs-Lewy (CFL) condition to be satisfied, that is The left-hand side in Eq. ( 29) is referred to as the Courant number (Courant et al., 1967).The domain of dependence of a point {x 0 , t 0 } is the set of points {x, t} upon which the solution u i (x 0 , t 0 ) has a dependency.The CFL condition ensures that the domain of dependence of the exact solution to the partial differential equation is contained within the domain of dependence of the numerical solution.If the CFL condition is not satisfied, then the exact solution depends on a larger set of points than the numerical solution.In such case, changing the state in certain points may change the exact solution but not the numerical solution, and therefore the numerical solution cannot converge to the exact one.

Validation of the method using benchmark problems
In the numerical simulations to follow, we use sets of values from Tab. 1 as the moduli of the phases.These values were taken as representative values within the characteristic range of the values that fits elastomers, see, e.g., Marckmann and Verron (2006), Carpi et al. (2008), Getz et al. (2017), and the references therein.In the numerical examples to follow, we set the Courant number to 0.9.We first test our scheme in Sec.4.1 against the analytical solution for the problem of a finite layer that is released from an initial shear strain, while bounded between two semi-infinite layers with different parameters.In Sec.4.2 we test it in a nonlinear problem of the scattering of an incident finite-amplitude shock wave from one half-space, impinging on an interface with a second half-space made of a different material.

Small-amplitude waves
Our first test case corresponds to a "fiber" of thickness l = 0.1 m with the property set 2 in Tab. ( 1), bounded by two semi-infinite "matrix" phases with the property set 1 in Tab. ( 1).The fiber is subjected to the infinitesimal initial strain field elsewhere. (30) This test case is similar to the test case by Bale et al. (2002), with the difference that we prescribe initial shear instead of compression, and use nonlinear phases.The solution is obtained using our finite volume method with a grid of 2000 cells per meter (∆X = 0.5mm).Fig. 3(a) shows the shear strain ǫ 2 as function of X 1 at t = 0 (dashed line) and t = 30 ms (solid line).As discussed by Bale et al. (2002), the initial strain is symmetrically separated into rightward and leftward propagating waves that repeatedly reflect and refract at the interfaces, thereby generating a train of (in this case shear) waves in the semi-infinite media.Note that ǫ 1 remains zero everywhere as it should in this uncoupled linear limit.The resultant waves are compared with the analytical solution for their length and velocity, see, e.g., the derivation by Shmuel and Moiseyev (2020).We first present in panel (b) the location of each strain local peak in the domain X 1 > 0.05 as function of time.We observe that lines are linear with the constant slop 9.98 ms −1 , which implies a 0.2% error with respect to the analytical value in the limit of small strains, namely, L = 10 ms −1 .Next, we compare the wavelengths in the numerical simulation to the analytical prediction.According to Eq. ( 46) of Shmuel and Moiseyev (2020), the allowed wavelength is λ by the peak-to-peak distance in our numerical simulation is exactly 0.1 m between every adjacent peaks.The distance between the first and second peaks is illustrated for example in Fig. 3(a).

Finite shock scattering at the interface between two half-spaces
We consider the case of reflection and transmission of an incident shock wave of finite amplitude at an interface between two elastic half-spaces.This nonlinear benchmark problem is chosen since we are able to obtain for it numerical solutions using standard root-finding algorithms-and specifically Newton's method-which will be compared with our finite-volume method.We use sets 1 and 2 for the moduli of phase m (left half-space) and phase f (right half-space, respectively.The left half-space is subjected to a shock wave for which the pre-shock state is unstrained and at rest, and its post-shock state is (ǫ 1 , ǫ 2 , v 1 , v 2 ) ≈ (−0.2, 2.3, 2.6, −32.4).These initial conditions were chosen to generate four shock waves, for which we obtain numerical solutions using Newton's method, as described by Ziv and Shmuel (2019) and in Appendix A. Fi g. 4(a) shows the distribution of ǫ 1 (upper panel) and ǫ 2 (lower panel) across −0.4 m < X 1 < 0.6 m at t = 10 ms.The black line corresponds to the numerical solution of the exact equations.To test the convergence of our scheme, we compare solutions using three different grids, namely, ∆X = 20 (orange triangles), 5 (red diamonds) and 2.5 mm (blue circles).It is clear that upon refining the grid, our scheme converges to the "exact" solution.In order to quantify the rate of convergence, we calculate the 1-norm of the error E between the two solutions, defined by where q is the "exact" solution.Fig. 4(b) shows the 1-norm of the errors of ǫ 1 (circles) and ǫ 2 (diamonds) as function of the cell length ∆X for each one of the three grids depicted in Fig. 4(a).The slope of the norms are around 1 in a logarithmic scale, indicating that the rate of convergence of the numerical solution is approximately of first order.This is expected since the second order correction of the method is neglected in the vicinity of discontinuities, where this solution consists of discontinuities only.

Numerical experiments of coupled waves in nonlinear laminates
We consider an infinite laminate composed of two alternating m and f compressible Gent layers with an equal length of H (m) = H ( f ) = 1 cm.We use sets 3 and 4 in Tab. 1 for the moduli of phases m and f , respectively.The laminate is subjected to the initial strain field where ǫ (I) i , A i and w are prescribed quantities, and the values of i depend on the type of displacements considered, as described next.Coupled axial and transverse displacements.-Here,i takes the values 1 and 2, and we set Fig. 5 shows ǫ 1 and ǫ 2 as functions of X 1 , using a grid of 2000 cells per meter (∆X = 0.5mm), such that each layer is discretized to 20 cells.Cyan, red, green and blue lines correspond to t = 5, 30, 55 and 80 ms, respectively.Remarkably, the rightward propagating wave5 is separated into a train of vector solitary waves, maintaining both their profile of the axial and shear components at the different periodic cells.This observation of vector solitary waves in Gent laminates follows the observation of scalar pressure solitary waves in nonlinear laminates by LeVeque (2002b) and LeVeque and Yong (2003b).Interestingly, we observe that the width of each vector solitary waves is approximately ten layers, similarly to the width reported by LeVeque and Yong (2003b) in the scalar case.While not shown here, we note that the number of generated solitary waves increases for greater values of w, as is the case for the KdV solitary waves.The dependency of the vector solitary wave velocity on the strain amplitude is studied in panel (b).Specifically, we show the velocity against the maximal value of ǫ 1 in the middle of phase m over the time that the solitary wave traverses a unit cell.We observe that the velocity is a monotonically increasing function of the strain amplitude.This dependency is in opposite of the dependency of the vector solitary waves discovered by Deng et al. (2017) in the discrete mechanical system they conceived, and similar to the dependency of the KdV solitons and the dependency observed by 5Owing to the symmetry of the problem, a mirrored wave that is propagating to the left is also generated.

LeVeque and Yong (2003b).
Coupled displacements in two dimensions.-Inthis case i = 2 and 3, and we set Fig. 6 shows ǫ T and θ as functions of X 1 , using a grid of 2000 cells per meter.Cyan, red, green and blue lines correspond to t = 20, 90, 160 and 230 ms, respectively.We observe that the transformation (18) proves useful also here, as it exposes two uncoupled polarizations of different nature.Specifically, its shows a linear polarization associated with ǫ T that breaks up into a train of three solitary waves with similar shape and different magnitude and velocity.The distance between each peak and its following peak is increasing in time, implying that the velocity is higher at higher strains.This is quantified in panel (b), where the velocity of the linearly polarized solitary waves is plotted against the maximal value of ǫ T at the middle of phase m over the time that the solitary wave traverses a unit cell.The second polarization associated with θ which lags behind the solitary waves is circular, and propagates effectively as a linear wave as it does not break up into θ-dependent waves; while not shown here for brevity, this occurs independently of the width of the initial localized strain w.Since the circular polarization is of linear waves, and it propagates independently of the linearly polarized solitary waves, the latter are scalar solitary waves and not vector solitary waves, even though they propagate through two coupled components of the displacement vector field.

Summary
We have developed a designated scheme to numerically solve the equations that govern elastic waves with two coupled components of finite amplitude in laminates made of nonlinear layers.Our scheme is based on two main elements.The first one is loaned from LeVeque (1997), whose finite-volume method utilizes the solution of a Riemann problem at the interface between grid cells to solve nonlinear hyperbolic systems that are not in conservation form.The second element is the flux-based wave decomposition of LeVeque (2002b) and Bale et al. (2002).Our extension required accounting for the generation of additional waves with respect to the acoustic problems of single stress component that were addressed by LeVeque (2002b) and Bale et al. (2002).This was carried out using a suitable matrix formulation which captures the coupling between the different components the displacement and stress fields.We have specifically addressed two cases, namely, a motion with coupling between its axial and transverse components, and a motion with two coupled transverse components.
We first tested our method using two benchmark problems.The first problem is linear, as there the initial strain is infinitesimal, and indeed our scheme recovered the analytical solution for the velocity and length of the generated waves.The second test case is of nonlinear waves and therefore more challenging.Specifically, we considered an incident shock of finite amplitude that strikes an interface between two half-spaces made of different nonlinear materials.The constitutive response of the half-spaces is described by the compressible Gent model (Gent, 1996, Lopez-Pamies andCastañeda, 2007); in our our previous work (Ziv and Shmuel, 2019) we have demonstrated that this model is capable of capturing shear shock waves and tensile-induced shocks phenomena in soft materials, which were observed experimentally by Catheline et al. (2003), Espíndola et al. (2017), andNiemczura andRavi-Chandar (2011), respectively.For this problem, numerical solutions using standard root-finding algorithms are accessible, and a comparison between such solutions using Newton's method and our method has shown an excellent agreement between the two solutions.This was the first numerical experiment of shock scattering between two half-spaces in finite elastodynamics with two coupled components.
Subsequently, we have applied our scheme in a numerical experiment of finite-amplitude waves with two coupled components in an initially strained periodic laminate made of two alternating compressible Gent layers.In the case of coupled axial and transverse displacements, our experiment revealed the generation of vector solitary waves.To the best of our knowledge, this is the first observation of vector solitary waves within the framework of continuum solid mechanics.Our observation was preceded by the first construction of vector solitary waves in discrete mechanical systems by Deng et al. (2017Deng et al. ( , 2018Deng et al. ( , 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 vectorial mechanical waves are slower at higher amplitudes, the vector waves in our continuum laminated model are faster at higher amplitudes.Therefore, the vector solitary waves here are more similar to KdV solitons and the acoustic solitary waves analyzed by LeVeque and Yong (2003b).
In the case of a coupling between two displacement components in the plane of the layers, our numerical experiment revealed the generation of a linearly degenerate wave of circular polarization that lags behind a train of solitary waves with linear polarization.Here again, the solitary waves are faster at higher amplitudes.
This paper establishes a starting point for several future works, among which we list the development of related higher-other methods such as the weighted essentially non-oscillatory (WENO) finite difference, see the survey by Shu (2016) and the references therein; methods and studies concerning higher-dimensional problems (LeVeque, 1997, Bale et al., 2002, Quezada de Luna and Ketcheson, 2014, Berjamin et al., 2019); analytical investigations using homogenization approaches (Andrianov et al., 2011, LeVeque and Yong, 2003b, Quezada de Luna and Ketcheson, 2014); and the introduction of a kinematic split between the isochoric and volumetric parts of the motion into the method, which is useful when the volumetric stiffness is several orders of magnitude greater than the shear stiffness (Holzapfel, 2000, Peyraut et al., 2007).

Figure 2 :
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 i and i + 1.

Figure 3 :Figure 4 :
Figure 3: (a) The distribution of ǫ 2 as function of X 1 at t = 0 (dash) and t = 30 ms (solid).(b) The location of each strain local peak in the domain X 1 > 0.05 as function of time.(a)(b)

Figure 5 :
Figure 5: (a) The strain distribution in the studied Gent laminate when subjected to axial-transverse strain.(b) The vector solitary wave velocity versus the maximal value of ǫ 1 in the middle of phase m over the time that the wave traverses a unit cell.

Figure 6 :
Figure 6: (a) The distribution of ǫ T and θ in the studied laminate when subjected to transverse-transverse strain.(b) The linearly polarized solitary wave velocity versus the maximal value of ǫ T at the middle of phase m over the time that the wave traverses a unit cell.

Table 1 :
Sets of Gent moduli used in the numerical simulations, chosen from the characteristic range of the values that fits elastomers.