A method for computation of discontinuous wave propagation in heterogeneous solids: basic algorithm description and application to one‐dimensional problems

An explicit integration algorithm for computations of discontinuous wave propagation in heterogeneous solids is presented, which is aimed at minimizing spurious oscillations when the wave fronts pass through several zones of different wave speeds. The essence of the present method is a combination of two wave capturing characteristics: a new integration formula that is obtained by pushforward–pullback operations in time designed to filter post‐shock oscillations, and the central difference method that intrinsically filters front‐shock oscillations. It is shown that a judicious combination of these two characteristics substantially reduces both spurious front‐shock and post‐shock oscillations. The performance of the new method is demonstrated as applied to wave propagation through a uniform bar with varying courant numbers, then to heterogeneous bars. Copyright © 2012 John Wiley & Sons, Ltd.


INTRODUCTION
High-fidelity analysis of wave propagation problems in damping-free homogeneous solids is obtained by a combination of three ingredients: the discrete equations of motion with regular uniform and equal mesh lengths, an explicit integrator with no numerical damping and minimum phase error, and the integration step size that satisfies the Courant number (C r D ct=x) to be unity, where .c, t, x/ are the speed of the sound of the material, the time step and the characteristic element length (or mesh size), respectively. In practice, although this ideal combination is difficult to realize, an acceptable compromise can often be achieved for discontinuous wave propagation analysis of homogeneous solids. Recently, a series of new synthesized materials including composites and metamaterials can be categorized as heterogeneous materials and are primary candidate materials in the design of optical and acoustic sensors, transportation systems and industrial machinery equipment. It is generally recognized that wave propagation analysis tools primarily developed for homogeneous solids become inadequate for heterogeneous materials, so much so that we are back to rely heavily on experiments, instead of complementing with simulation to a great extent, 623 for material characterization of most of new synthesized materials for dynamic applications. It is with this background that we are motivated to develop an algorithm for wave propagation through heterogeneous materials.
For the successful tracing of wave propagations through heterogeneous media, three approaches may be explored: (1) Employ a uniform step size for all the heterogeneous zones and apply a compensation strategy that would minimize spurious oscillations for zones whose effective Courant number fall into C r D ct=x << 1.
(2) Design an effective spatial and/or temporal filter that alleviates spurious spillovers whenever waves are passing through and reflecting the heterogeneous interfaces. An example of this approach involves a decimation of spurious reflecting waves by constructing a set of windowed Fourier filters. Although proven useful for linear homogeneous problems, such approaches are of limited use for nonlinear wave propagation problems. (3) Employ variable step sizes for each of heterogeneous zones to minimize spurious oscillations inherent for all the existing integration schemes when the step size is different from the ideal one, namely C r D ct=x D 1. This necessitates an adoption of subcycling integration strategies with considerable implementation complexities.
For clarity, we introduce the following definitions: Post-shock oscillations refer to the spurious oscillations that occur behind the shock fronts. We will also use synonymously after-shock oscillations.
Front-shock oscillations refer to the spurious oscillations that occur in front of the discontinuities. We will also use synonymously before-shock oscillations.
Common to all of these three approaches is the fact that, to a varying degree, all the existing explicit numerical schemes suffer from engendering post-shock spurious oscillations unless one integrates each zone with the ideal stepsize, C r D ct=x D 1. When an implicit integration scheme is used, both post and front shocks may occur, depending on the nature of spatial discretizations. By post-shock oscillations, it is meant to refer to deleterious oscillations trailing the stress waves and/or velocity discontinuities. It is well known that post-shock oscillations are triggered by the mismatch of the physics of wave propagation speed across the adjacent two discrete nodes and  the inherent numerical coupling stemming from the discretization of the field equation. It is observed that regardless of the step size, the wave front moves one node per step. It so happens that, with the ideal step size, C r D ct=x D 1, the wave front moves exactly one node per step. For example, when one takes the step size of C r D ct=x D 1=2, even though it should physically take two steps for the wave front to move one node, the wave front advances two nodes. This is illustrated in Figure 1, where a unit displacement is introduced at the bar center and let it propagate. Observe that at time t D t, the wave fronts should be at nodes 5 and 7 as shown in the right part of Figure 1.
On the other hand, when integrated with the step size t D 0.5x=c, the wave fronts are already at nodes 4 and 8 although the magnitude is relatively small. At time t D 3t , the correct wave fronts should be at nodes 3 and 9 as shown in the right part of Figure 1. However, when integrated with the stepsize t D 0.5x=c as shown in the left part of Figure 1, the wave fronts are erroneously at nodes 2 and 10, thus creating numerical dispersion. This is explained later. Consider the well-known three-point discrete approximation of @ 2 w @x 2 (or linear two-node approximation in the finite element method) in one dimensional wave equation in an elastic bar as shown: Continuum wave equation : @ 2 w @t 2 c 2 @ 2 w @x 2 D f .x, t/, 0 6 x 6 L, 0 6 t 6 t max with initial/boundary conditions: Semi-discrete equation at j th node: Semi-discrete equation at .j C 1/th node: is approximated by the spatial central difference scheme; w.x, t/ is the displacement; .x, t/ are material point in the bar and time, respectively; c is the speed of the sound of the material of the elastic bar; and, .w l .t /, w r .t / are the boundary displacements at the left and right ends of the bar, respectively.
The preceding semi-discrete equations reveal that if the wave front is at the j th node at time t , that is, w j 6 D 0, then it will give non-zero acceleration at the .j C 1/-node as can be seen from the second discrete equation for .j C 1/-node regardless of the step size. In other words, the physical shock front is only halfway between the next node when C r D ct=x D 1=2. Nonetheless, the numerical scheme gives a false pretense that the wave front arrives at the next node. It is this disparity between the actual physics and the numerical artifice that is responsible for post-shock spurious oscillations as illustrated in Figure 2.
The preceding observations have led us to ponder the following conjecture: Can one develop an explicit algorithm that captures discontinuous waves with front-shock oscillations while minimizing post-shock oscillations? If this is possible, then perhaps one may be able to alleviate both the front and post-shock oscillations by a judicious averaging of the front-shock inducing and post-shock inducing algorithms. We offer a preliminary indication that the answer to the preceding question is affirmative as shown in Figures 3 and 4.
The field of shock capturing algorithms in solids and fluids abounds with a rich body of literature proposed by many researchers over the years. Beginning with the seminal works of Courant, Friedrichs and Lewy [1], Lax [2], Lax and Wendroff [3], and the classical text of Rightmyer and Morton [4], shock capturing algorithms for wave propagation have advanced in several avenues: the finite difference method [5][6][7][8][9][10][11], front tracking algorithms [12][13][14][15][16], finite element spatial discretization with the finite difference in time [17][18][19][20][21][22][23][24][25][26], space-time treatment methods [27][28][29]   Galerkin's methods [30][31][32][33], oscillations filtering by postprocessing [34,35], and recently, variational construction [36][37][38][39][40]. As for the time discretization, for wave propagation in solids, the central difference method is most widely used for explicit treatment in time and the Newmark average acceleration method for implicit time integration [41]. Interested readers are referred to a special edition of finite element-based computational methods in Wave Motion [21] and Zuazua [10] for finite difference-based methods up to the mid-2000s, among others. In the context of the literature cited earlier, the present proposed method may be viewed as a Lagrangian front tracking method that is tailored for wave propagation in solids. Specifically, the present method employs the finite element-based spatial discretization that is analogous to centered spatial finite difference discretization. The time discretization is a combination of the central difference method and a pushforward-pullback interpolation scheme for wave-front tracking. It will be shown that the proposed method significantly alleviates spurious oscillations when the integration step sizes are different from the critical stepsize, namely ct=x < 1. Unlike methods that filter spurious oscillations by post processing (see, e.g., [34,35]), the present method does not require filtering by post processing, which is essential for applications to nonlinear problems.  The classical central difference method is well-known to trigger post-shock spurious oscillations, while having a built-in ability of filtering the front-shock oscillations. On the other hand, the proposed pushforward-pullback front tracking algorithm engenders front-shock oscillations while filtering out the post-shock oscillations. These complementary characteristic features, when combined together, significantly alleviate both spurious post-shock and font-shock oscillations. The rest of the paper is organized as follows. Section 2 introduces the basic concept of the present algorithm. The present algorithm, first, employs a pushforward explicit integrator that pushes the wave fronts, regardless of the stepsize taken (t) relative to the critical stepsize (t c ), to the next wave front nodes as if the stepsize taken is t c . A pullback integrator is then employed to adjust the wave fronts to be where they should be when˛D t=t c < 1. The net effect of the pushforward and pullback integration is demonstrated as applied to a box wave propagation. Also shown is the intrinsic feature of the central difference method that triggers post-shock spurious oscillations. Section 3 introduces the present method that combines the pushforward-pullback integrator and the central difference method and applies the new method to discontinuous wave propagation through a homogeneous bar. The development of the present algorithm applicable to heterogeneous solids is presented in Section 4, which involves a careful treatment of interface nodes. Section 5 offers stability and accuracy assessments of the present algorithm. The analysis results indicate that the new algorithm recovers the identical stability limit of that of the central difference method when t D t c , and maintains sufficient stability for reduced stepsizes, namely˛< 1. The present algorithm does induce numerical damping, which remains small, 0.7% when ! max t D 1. As for frequency errors, the present method does increase the phase error about twice that of the central difference method when˛D 0.5 when no dispersion minimization is introduced. Section 6 reports numerical performance of the present method as applied to heterogeneous bars, with substantially reduced spurious oscillations compared with existing shock capturing algorithms such as the central difference method, the fourth-order Runge-Kutta method and their allied methods. Finally, a summary of the present study is offered in Section 7.

CONSTRUCTION OF AN ALGORITHM FOR INDUCING FRONT-SHOCK OSCILLATIONS
We have shown in the Introduction that wave fronts advance one discrete node per time step no matter how small the step size is relative to the ideal time step size, namely C r D ct c =x D 1.

627
It is recalled that when the step size is less than the ideal one, t < x=c, the coupling term for the j th equation, forces part of the wave energy to be allocated at the .j C 1/-node. This happens in spite of the fact that the wave front is impossible to reach from the j th node to the .j C 1/-node within one step whenever t < x=c. Thus, the coupling terms in effect play the role of artificially stretching the wave speed to where t c corresponds to the ideal step size that makes the Courant number to be unity. It is emphasized that this artificially stretched wave speed (c artificial ) gives rise to the post-shock spurious oscillations. This observation suggests that if one artificially increases the step size that is larger than the actual step size, then perhaps the spurious oscillations may exhibit in front of the discontinuities instead of behind the discontinuities. To this end, we explore the following combination of extrapolation and interpolation.
First, we discretize the one-dimensional continuum bar equation (1) by the standard finite element spatial discretization procedure to result in the following semi-discrete equations of motion as where .M, K/ are the mass and stiffness matrices; .w, f/ are the displacement and applied force vectors; and the superscript dot designates time differentiation. Second, assuming that the displacement (w) and the velocity ( P w) at time t D nt are available, we obtain by an extrapolation (or by an explicit integration formula) the displacement at time t nCc D t n C t c as shown in Figure 5.
Here, for simplicity, we have employed the standard central difference method. Third, we compute the acceleration at the extrapolated step .n C c/: Fourth, we interpolate the displacement at .nC1/ by utilizing the displacements and accelerations at the nth step and the extrapolated (n C c/th step: which may be viewed as a discrete pullback operation. Once the displacements, w nC1 , are updated, the velocity vector ( P w nC1 ) is updated via Figure 5. Extrapolation of the displacement w.x, t/ at step .nCc/ followed by interpolation at the next step .n C 1/. Here, we employ the shock front at .n C c/ for the critical imagined stepsize t c . Then, we pull back the wave front to .n C 1/ that is commensurate with the actual stepsize t.

Remark 1
Note that, when˛D 1.t D t c /, we havě Hence, when one integrates with the step size of t D t c D x=c, the present front-shock inducing algorithm (6) coincides with the central difference method. For this special case, stress waves and/or discontinuous forcing functions do not trigger spurious oscillations.

Remark 2
The algorithm presented in (4)- (7) requires two function evaluations to compute R w nCc (5) and R w nC1 (7). Hence, one is tempted to obtain R w nC1 by interpolating R w n and R w nCc : to reduce function evaluations per step from two to one. As illustrated in Figure 6, (7) induces front-shock oscillations whereas (9) does not.

DEVELOPMENT OF PRESENT SHOCK CAPTURING ALGORITHM IN SOLIDS
It was alluded to in the Introduction that almost all of the existing explicit algorithms trigger post-shock oscillations. To this end, we illustrate box-wave propagation along a one-dimensional bar obtained by two widely used integrators: the explicit fourth-order Runge-Kutta method and the central difference method when t c < x=c as shown in Figure 7. It should be noted that the trapezoidal rule (known also as the mid-point implicit method and the Newmark algorithm ( D 0.5,ˇD 0.25)) is shown to be incapable to treat momentum discontinuity in [36]. Hence, we would not explore the usage of implicit methods any further.  (7) and (9) for inducing front-shock oscillations with D 0.5. With one function evaluation, formula (9) fails to induce front-shock oscillations. Two function evaluations per step given by formula (7) successfully induces front-shock oscillations.  We now focus on the exploitation of two explicit methods: the front-shock inducing method described in the preceding section as described in (4)-(7) and a version of the standard central difference method given by For convenience, we summarize the front-shock inducing method presented in the preceding section as follows: Let us now combine the two methods (10) and (11) to obtain Figure 8 shows the box-wave tracing by the new method (12), which, when compared with either that by the post-shock oscillating central difference method (10) alone or by the front-shock inducing method (11) alone as shown in Figure 9, exhibits a dramatic improvement. Although not reported herein, we have performed additional numerical experiments for all ranges of the Courant numbers, namely 0 < ct=x 6 1.0, which proved that the new method consistently improves shock capturing ability compared with existing explicit methods. Encouraged by the substantial improvements by the new method (12), we will evaluate the new method for wave propagation through heterogeneous bar. This is described in the next section.

PRESENT ALGORITHM FOR HANDLING HETEROGENEITY
For clarity, we consider two bars with significantly different wave speeds as shown in Figure 10. As the interface needs a special treatment, we present a step-by-step procedure utilizing the present algorithm (10)- (12).  Figure 10. Bar consisting of two heterogeneous materials. Figure 11. Extrapolation of the displacements .w 1 .x, t/, w 2 .x, t// at steps ..n C c1/, .n C c2// followed by interpolations at the next step .n C 1/ for heterogeneous problems.
First, we obtain the pushforward displacements for the two domains, w nCc1 1 , w nCc1 2 as follows (See Fig. 11) It should be noted that t c1 is the maximum pushforward timestep without triggering instability. Using these pushforward displacements, we obtain the accelerations R w nCc1 1 , R w nCc1 2 from the following partitioned equation set (see, e.g., [42] for details in the partitioned solution procedures): where (B 1 , B 2 ) are the Boolean matrices that extract the interface accelerations, (L f 1 and L f 2 ) are the extractor for the frame acceleration ( R w f ); ( 1 and 2 ) are the interface forces, respectively, as shown in Figure (12); and .w 1b and w 2b / are the interface boundary displacement for materials 1 and 2, respectively.
In the preceding equation set, the third and fourth equations state that the boundary displacements (w 1b , w 2b ) are related to the frame displacement (w f ). The fourth equation states that the interface force sum is equal to zero, stating Newton's third law.
It should be noted that computations of ( 1 , 2 , R w f ) require, for the one-dimensional bar, the solution of a .3 3/ matrix equation, and for multi-dimensional problems, a .3n g 3n g / matrix Heterogeneous bar parttioned w f w = 1b 1 w 1 w = 2b 2 w 2 t t Fixed end Figure 12. Heterogeneous bar partitioned. Here, we treat the displacement at the interface node, which is preserved as a frame displacement (w f ) and two distinct Lagrange multipliers are introduced to enforce the two independent interface compatibility conditions (cf., Equation 14). equation where n g is the size of interface displacement vector. The pushforward displacement and acceleration for the material (2) at the critical time step .n C c2/ is obtained by The pullback displacements at the time step .n C 1/ are thus obtained by In these equations, (t c1 , t c2 ) are the critical stepsizes for materials (1) and (2), respectively, with their corresponding wave speeds (c 1 , c 2 ). Second, we obtain the pullback accelerations: where p 1,2 are the momentum vector of two domains, L g is a Boolean operator that assembles the two domain displacements into the corresponding assembled displacement, and R w g is the assembled acceleration vector.
Third, we update the velocity for each domain via Fourth, we obtain the assembled velocity and displacement by where L C g is a pseudo-inverse of L g . Fifth, we compute the post-shock oscillating response via where the subscript g denotes the assembled quantities.

633
Sixth, use (12) to obtain the acceleration, velocity and the displacement, which we repeat for completeness: For partitioned computations of the two heterogeneous bars, we obtain the partitioned accelerations, velocity and displacement by This completes the computational sequence for a time step.

STABILITY AND ACCURACY ANALYSIS OF THE PRESENT ALGORITHM
It turns out that, for linear and linearized equations of motion, the stability of the present algorithm is dictated by the combination of the pushforward and pullback formulas for the stiff bar. This observation allows one to assess the stability of the present algorithm for heterogeneous cases by assuming that the entire bar consists of stiff materials. In other words, one can assess the stability by employing the present algorithm as specialized to a homogeneous bar (10)- (12). To this end, we obtain the following difference equation from (10)- (12): where ! k is the undamped frequency for the kth mode, T is the mass-normalized eigenvector. Focusing on the highest frequency of the system, the discrete characteristic equation of the mentioned difference equation is obtained by seeking w nC1 and w n in terms of w n 1 by w nC1 D 2 w n 1 , w n D 2 w n 1 (25) yielding for the highest frequency the following characteristic equation: For stability, the amplification factor must be bounded by which can be mapped onto the entire left-hand plane of the´-plane given by For stability: a > 0, b > 0, c > 0 via Routh-Hurwitz criterion (29) For˛6 1, we observe with D 1 2 and Â D 1 2 that the conditions of b > 0 and c > 0 are satisfied. This leaves to examine a: Figure13 shows the stability of limits of the present method plotted for stable zones versus the stepsize ratio (˛D t=t c ). In that figure, when˛D 1 (viz., t D t c ), the stability limit is given by When the stepsize () is reduced from the critical stepsize (t c ), the sufficient stable stepsize is given by This is shown as a straight line in Figure 13. Plotted on the same figure is the stability condition (30). Note that for all the values of stepsizes, the present method maintains stability with sufficient stability margins.
To assess the accuracy of the present method, the numerical damping and frequency errors are plotted in Figures 14 and 15 for the stepsize (˛D t=t c D 0.5). Whereas the central difference method (˛D 1) introduces no damping (j j D 1.0), the present method does introduce when < 1.0. Nevertheless, the amount of numerical damping is less than 7% when ! max t D 1.0, which is rather a large stepsize for tracing oscillatory signals.
Finally, the frequency error of the present method is plotted for the case of˛D 0.5 versus ! max t as shown in Figure 15. Note that the frequency is increased for the central difference method (˛D 1) whereas the present method with˛D 0.5 decreases the frequency. The absolute error of the present   method is about twice that of the central difference method. One may employ the mass averaging method that was successfully adopted in [36] to minimize dispersion errors, but was not pursued in the present study and left for subsequent investigations.

NUMERICAL EXPERIMENTS FOR WAVE PROPAGATION THROUGH HETEROGENEOUS BAR
To demonstrate the performance of the present heterogeneous algorithm as detailed in (13)-(22), we apply them to two example problems as discussed later. For comparison purposes, we have used the explicit variational integrator presented in [36], whose basic algorithm is an adaptation of the central difference algorithm with a feature for reducing dispersion errors.

Box wave applied at stiff, free-end and propagating into soft, fixed-end bar
The proposed explicit integration algorithm for handling heterogeneous bars (13)-(22) is applied to a bar consisting of a stiff bar on the left side and a bar made of soft material on the right side as shown in Figure 16. The ratio of wave speeds between the soft and stiff bar is 10, leading to the Courant number ratio to the same 10. A box wave is applied at the left end, and the velocity profiles along the bar versus time elapsed is plotted therein. The central difference method is compared with the variational explicit integrator presented in [36]. As can be seen, the central difference method exhibits considerable post-discontinuity oscillations. Although not reported in the present paper, the deleterious spurious oscillations engendered by the central difference method are observed for other cases. Hence, the performance of the central difference method is not reported in the case of heterogeneous bars. These results demonstrate that the proposed new explicit integrator yields robust wave tracing compared with the explicit variational integrator. Again, we have carried out comparisons, although not reported herein, with the fourth-order Runge-Kutta method in addition to the central difference method. Figure 17 compares the performance of the explicit variational integrator [36] and the present method (13)- (22). Notice that at time t D 85, the wave front is still confined within the stiff bar on the left side as the interface of the two bars is at 100 on the bar coordinate. Here the performance of the two methods are similar.
When the wave front moves into the soft bar as shown in Figure 18, the explicit variational integrator exhibits appreciable spurious oscillations. The spurious oscillations do not subside away as the wave front moves further into the soft bar as evidenced in Figure 19.

Box wave applied at soft, free-end and propagating into stiff, fixed-end bar
This time, we place the soft bar on the left end and subject to the same box wave applied at the left free end. When the step size corresponding to the Courant number for the stiff portion of the bar to be C r D c stiff t=x D 0.85 is chosen, the corresponding Courant number for the soft bat becomes C r D 0.085, far below the desirable value of unity. The significant spurious oscillations engendered by the variational explicit integrator [36] is evident as shown in Figure 20. Figure 21 reports the velocity profiles at later times, t D 1105 and t D 1700, corresponding to the wave front passing through the interface and reflecting from the fixed end, respectively. The results indicate that the new method mitigates the spurious oscillations for waves propagating through heterogeneous bars. It is noted that the central difference method, although not reported in Figures 20  and 21, yields far worse spurious oscillations than the variational integrator [36].

Box wave applied at soft-stiff-stiffer and stiffer-stiff-soft heterogeneous bars
We have applied to a heterogenous bar consisting three differing materials and/or differing mesh sizes. The geometry of the bar consisting of three differing materials are shown in Figure 22. Two cases are examined: soft-stiff-stiffer materials from the left end to the right end, and stiffer-stiff-soft materials from the left end to the right end. A box wave is applied whose durations correspond to cover half of the first bar for both cases. Specifically, we have applied a unit applied force at the left end in Figure 22 given by where .L left , c left / are the bar length and the speed of propagation of the left end bar; and, .T , t, N / are the time duration of applied constant force, t is the integration step size, and N is the number of integration steps, respectively. In all of the four numerical experiments reported, we have used 100 elements for each of the three material regions, and the speed of propagation of these bar segments are shown in the subsequent figures. Figure 23 shows the response to the box wave applied to at the left end (the soft materials segment). The analytical wave front should be at the middle of the second bar segment, denoted as stiff materials in Figure 23). Because of inherent dispersion effects of numerical integrators when the step size is less than the ideal one, namely C r < 1.0, the numerical wave front is approximate at three quarter of the second bar segment. Note the excessive spurious oscillations for response obtained by the central difference method compared with the response obtained by the present method. This is because the step size is determined by the stiffer portion of the bar with C r stiffer D 0.8, which makes that of the soft segment to be C r soft D 0.080, a hundredth of the stiffer segment, causing excessive spurious oscillations. Figure 24 shows the wave pattern when the analytical wave front should be at the middle of the stiffer bar segment. Once again, the numerical dispersion gives rise to the numerical wave front near at the fixed end. Observe that the present method continues to perform in a robust manner compared with the central difference method. In particular, the refracting waves back into the soft bar segment cause significant oscillations with the central difference method. In addition, the response of the middle segment of the bar consists of both refracting waves and as well as the advancing waves. One notices a smooth yet spurious oscillations by the central difference method.
In Figure 25, the bar is arranged in an opposite manner: stiffer-to-stiff-to-soft arrangement. This is indicated in Figure 25. Note that when the analytical wave front is at the middle of the second bar segment (stiff bar portion), the central difference method experiences significant spurious oscillations both at the stiffer segment and stiff segment. Observe when the wave propagate from stiffer to less stiff (designated as stiff segment), there is far less numerical dispersion than for the case of a wave propagating from soft to the stiff materials as demonstrated in Figures 23 and 24.  Finally, Figure 26 shows the propagating wave pattern consisting of both advancing and refracting waves when the analytical wave front is at the middle of the soft bar segment.
Once again, the present method offers a robust solution compared with the central difference method. Note that there is very little wave dispersion on the wave front as indicated in Figure 26. This and a similar phenomenon observed in Figure 25 confirms that, when a wave propagates from the stiff materials into soft materials, one expects very little numerical dispersion, whereas when waves are propagating from soft into stiff materials, significant (hence so far unavoidable) dispersions may be exhibited.

DISCUSSIONS
A method for tracing wave propagation in solids is presented, which obviates spurious oscillations when the integration step sizes are different from the critical stepsize, namely ct=x < 1. The present method consists of a combination of two explicit methods: the classical central difference  method that can trigger post-shock spurious oscillations, and a new pushforward-pullback front tracking algorithm that engenders front-shock oscillations. A closer observation of these two methods reveals that the central difference method has an ability to filter out the front-shock oscillations, whereas the new pushforward-pullback front tracking method filters out the post-shock oscillations. A weighted combination of the two methods has led to the development of the present method, initially for the homogeneous materials (10)- (12), and subsequently for heterogeneous materials (13)- (22). Hence, the present method can be used for computing wave propagation in heterogeneous solids for which one has to deal with a widely varying Courant numbers for a given integration step size. For clarity of presenting the essence of the present method, we have not included additional embellishments such as mass averaging to improve dispersion errors and numerical damping to chop off blips before and after the discontinuities (see Figures 4 and 9). We have worked on these embellishments in the development of variational integrators [36], which we intend to introduce in solving two and three-dimensional heterogeneous wave propagation problems.
Whereas the proposed method has shown to be effective for mitigating spurious oscillations when C r D ct=x < 1 for one-dimensional problems, we submit that the present paper is an initial report of our studies for the development of robust shock capturing algorithms for multi-dimensional heterogeneous solids. For example, we have not dealt with minimizing dispersion errors as well as smoothing the shock profiles. We plan to report on these aspects in conjunction with an adaptation of the present method for shock propagation in multi-dimensional and multi-scale problems.
A potential adaptation of the present method is for discontinuous wave propagation in fluid and gas dynamics problems. Another is to apply it to subcycling integration of the dynamics responses of heterogeneous materials. As the required extension needs further refinements and additional considerations, these and additional related issues are being actively investigated.
Finally, a formalization of the present method cast either in variational integrator or fractional step forms would be desirable so that one can exploit symmetry and symplecticity properties [43,44] for long-term wave propagation calculations. These and other issues remain open problems allied with the present method.