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

K.C. Park, S.J. Lim, H. Huh

## ► To cite this version:

K.C. Park, S.J. Lim, H. Huh. A method for computation of discontinuous wave propagation in heterogeneous solids: basic algorithm description and application to one-dimensional problems. International Journal for Numerical Methods in Engineering, 2012, 91, pp.622 - 643. 10.1002/nme.4285. hal-01443190

HAL Id: hal-01443190

https://hal.science/hal-01443190

Submitted on 22 Jan 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A method for computation of discontinuous wave propagation in heterogeneous solids: basic algorithm description and application to one-dimensional problems

K. C. Park<fnref n="1" />,<fnref n="2" />, S. J Lim<fnref n="3" /> and H. Huh<fnref n="3" />

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.

**KEY WORDS:** explicit integrator, minimal spurious oscillations, heterogeneous bar, propagation of stress waves

## 1. 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 ($Cr = c\Delta t/\Delta x$) to be unity, where ($c, \Delta t, \Delta 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,

<footnote n="1">¹ Department of Aerospace Engineering Sciences, University of Colorado at Boulder, Boulder, CO 80309-429, USA</footnote>

<footnote n="2">² Division of Ocean Systems Engineering, KAIST, Daejeon 305-701, Korea</footnote>

<footnote n="3">³ Division of Mechanical Engineering, School of Mechanical, Aerospace and Ocean Engineering, KAIST, Daejeon, 305-701, Korea</footnote>
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 $Cr = c\Delta t/\Delta x \ll 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 $Cr = c\Delta t/\Delta x = 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, $Cr = c\Delta t/\Delta x = 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

<figure><img src="image_2.png" /><figcaption>Figure 1. Wave propagation with dispersion caused by disparity of wave speed and the time step size with $Cr = c\Delta t/\Delta x = 1/2$ (left) and without dispersion $Cr = c\Delta t/\Delta x = 1$ (right). The central difference method is used to trace the wave fronts.</figcaption></figure>
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, $Cr = c\Delta t/\Delta x = 1$, the wave front moves exactly one node per step. For example, when one takes the step size of $Cr = c\Delta t/\Delta x = 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 = \Delta 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 $\Delta t = 0.5\Delta x/c$, the wave fronts are already at nodes 4 and 8 although the magnitude is relatively small. At time $t = 3\Delta t$, 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 $\Delta t = 0.5\Delta x/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 $\frac{\partial^2 w}{\partial x^2}$ (or linear two-node approximation in the finite element method) in one dimensional wave equation in an elastic bar as shown:

$$ \text{Continuum wave equation: } \frac{\partial^2 w}{\partial t^2} - c^2 \frac{\partial^2 w}{\partial x^2} = f(x,t), \quad 0 \le x \le L, \quad 0 \le t \le t_{\max} $$

$$ \text{with initial/boundary conditions: } \frac{\partial w(x, 0)}{\partial t} = \dot{w}(x, 0), \quad 0 \le x \le L \\ w(0, t) = w_l(t), \quad w(L, t) = w_r(t) \\ \Downarrow $$

$$ \text{Semi-discrete equation at } j\text{th node: } \frac{\partial^2 w_j}{\partial t^2} - c^2 \frac{w_{j-1} - 2w_j + w_{j+1}}{\Delta x_j} = 0 $$

$$ \text{Semi-discrete equation at } (j+1)\text{th node: } \frac{\partial^2 w_{j+1}}{\partial t^2} - c^2 \frac{w_j - 2w_{j+1} + w_{j+2}}{\Delta x_{j+1}} = 0 \quad (1) $$

where $(\frac{\partial^2 w}{\partial x^2})$ 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 \neq 0$, then it will give non-zero acceleration at the $(j+1)$-node as can be seen from the second discrete equation for $(j+1)$-node regardless of the step size. In other words, the physical shock front is only halfway between the next node when $Cr = c\Delta t/\Delta x = 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–11], front tracking algorithms [12–16], finite element spatial discretization with the finite difference in time [17–26], space-time treatment methods [27–29],
<figure><img src="image_3.png" /><figcaption>Figure 2. Post-shock spurious oscillations of a box wave passing through a homogeneous bar with step size of $Cr = c\Delta t/\Delta x = 1/2$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Front-shock spurious oscillations of a box wave passing through a homogeneous bar with step size of $Cr = c\Delta t/\Delta x = 1/2$.</figcaption></figure>

Galerkin's methods [30–33], oscillations filtering by postprocessing [34, 35], and recently, variational construction [36–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 $c\Delta t/\Delta 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.
<figure><img src="image_5.png" /><figcaption>Figure 4. Solution of a front-shock inducing algorithm and the post-shock oscillation algorithm with step-size of $Cr = c\Delta t/\Delta x = 1/2$, obtained by averaging the results of a front-shock inducing algorithm and a post-shock inducing algorithm at each time step.</figcaption></figure>

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 front-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 ($\Delta t$) relative to the critical stepsize ($\Delta t_c$), to the next wave front nodes as if the stepsize taken is $\Delta t_c$. A pullback integrator is then employed to adjust the wave fronts to be where they should be when $\alpha = \Delta t/\Delta 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 $\Delta t = \Delta t_c$, and maintains sufficient stability for reduced stepsizes, namely $\alpha < 1$. The present algorithm does induce numerical damping, which remains small, 0.7% when $\omega_{max}\Delta t = 1$. As for frequency errors, the present method does increase the phase error about twice that of the central difference method when $\alpha = 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.

## 2. 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 $Cr = c\Delta t_c/\Delta x = 1$.
It is recalled that when the step size is less than the ideal one, $\Delta t < \Delta x/c$, the coupling term for the jth equation, forces part of the wave energy to be allocated at the (j + 1)-node. This happens in spite of the fact that the wave front is impossible to reach from the jth node to the (j + 1)-node within one step whenever $\Delta t < \Delta x/c$. Thus, the coupling terms in effect play the role of artificially stretching the wave speed to

$$c_{\text{artificial}} = c \frac{\Delta t_c}{\Delta t}, \quad \Delta t_c = \Delta x/c \qquad (2)$$

where $\Delta 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_{\text{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

$$M\ddot{\mathbf{w}}(t) + K\mathbf{w}(t) = \mathbf{f}(t), \quad \mathbf{w}(0) = \mathbf{w}^0, \quad \dot{\mathbf{w}}(0) = \dot{\mathbf{w}}^0 \qquad (3)$$

where $(M, K)$ are the mass and stiffness matrices; $(\mathbf{w}, \mathbf{f})$ are the displacement and applied force vectors; and the superscript dot designates time differentiation.

Second, assuming that the displacement $(\mathbf{w})$ and the velocity $(\dot{\mathbf{w}})$ at time $t = n\Delta t$ are available, we obtain by an extrapolation (or by an explicit integration formula) the displacement at time $t_{n+c} = t_n + \Delta t_c$ as shown in Figure 5.

$$\mathbf{w}^{n+c} = \mathbf{w}^n + \Delta t_c \dot{\mathbf{w}}^n + \frac{\Delta t_c^2}{2} \ddot{\mathbf{w}}^n \qquad (4)$$

Here, for simplicity, we have employed the standard central difference method.

Third, we compute the acceleration at the extrapolated step $(n + c)$:

$$\ddot{\mathbf{w}}^{n+c} = \mathbf{M}^{-1}(\mathbf{f}^{n+c} - K\mathbf{w}^{n+c}) \qquad (5)$$

Fourth, we interpolate the displacement at $(n+1)$ by utilizing the displacements and accelerations at the nth step and the extrapolated $(n+c)$th step:

$$\mathbf{w}^{n+1} = \mathbf{w}^n + \Delta t_c \alpha \dot{\mathbf{w}}^n + \Delta t_c^2 \beta_1(\alpha) \ddot{\mathbf{w}}^n + \Delta t_c^2 \beta_2(\alpha) \ddot{\mathbf{w}}^{n+c} \\ \beta_1(\alpha) = \frac{1}{6}\alpha(1 + 3\alpha - \alpha^2), \quad \beta_2(\alpha) = \frac{1}{6}\alpha(\alpha^2 - 1), \quad \alpha = \frac{\Delta t}{\Delta t_c} \qquad (6)$$

which may be viewed as a discrete pullback operation.

Once the displacements, $\mathbf{w}^{n+1}$, are updated, the velocity vector $(\dot{\mathbf{w}}^{n+1})$ is updated via

$$\ddot{\mathbf{w}}^{n+1} = \mathbf{M}^{-1}(\mathbf{f}^{n+1} - K\mathbf{w}^{n+1}) \\ \Downarrow \\ \dot{\mathbf{w}}^{n+1} = \dot{\mathbf{w}}^n + \Delta t \gamma \ddot{\mathbf{w}}^n + \Delta t(1-\gamma)\ddot{\mathbf{w}}^{n+1}, \quad 0.5 \le \gamma < 1 \qquad (7)$$

<figure><img src="image_6.png" /><figcaption>Figure 5. Extrapolation of the displacement $w(x, t)$ at step $(n+c)$ followed by interpolation at the next step $(n+1)$. Here, we employ the shock front at $(n+c)$ for the critical imagined stepsize $\Delta t_c$. Then, we pull back the wave front to $(n+1)$ that is commensurate with the actual stepsize $\Delta t$.</figcaption></figure>
### Remark 1

Note that, when $\alpha = 1(\Delta t = \Delta t_c)$, we have

$$ \beta_1(\alpha) = \frac{1}{2}, \quad \beta_2(\alpha) = 0 \qquad (8) \\ \mathbf{w}^{n+c} = \mathbf{w}^{n+1} $$

Hence, when one integrates with the step size of $\Delta t = \Delta t_c = \Delta 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 $\ddot{\mathbf{w}}^{n+c}(5)$ and $\ddot{\mathbf{w}}^{n+1}(7)$. Hence, one is tempted to obtain $\ddot{\mathbf{w}}^{n+1}$ by interpolating $\ddot{\mathbf{w}}^n$ and $\ddot{\mathbf{w}}^{n+c}$:

$$ \ddot{\mathbf{w}}^{n+1} = (1 - \alpha)\ddot{\mathbf{w}}^n + \alpha\ddot{\mathbf{w}}^{n+c} \qquad (9) \\ \downarrow \\ \dot{\mathbf{w}}^{n+1} = \dot{\mathbf{w}}^n + \frac{1}{2}\Delta t (\ddot{\mathbf{w}}^n + \ddot{\mathbf{w}}^{n+1}) $$

to reduce function evaluations per step from two to one. As illustrated in Figure 6, (7) induces front-shock oscillations whereas (9) does not.

## 3. 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 $\Delta t_c < \Delta 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 ($\gamma = 0.5, \beta = 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.

<figure><img src="image_7.png" /><figcaption>Figure 6. Comparison of formulas (7) and (9) for inducing front-shock oscillations with $\gamma = 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.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7. Post-shock spurious oscillations by the Runge-Kutta fourth-order method and the central difference method when integrated with $\Delta t = 0.5\Delta x/c$.</figcaption></figure>

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

$$ 
\begin{aligned}
\mathbf{w}_{cd}^{n+1} &= \mathbf{w}^n + \Delta t \dot{\mathbf{w}}^n + \frac{\Delta t^2}{2} \ddot{\mathbf{w}}^n \\
\ddot{\mathbf{w}}_{cd}^{n+1} &= \mathbf{M}^{-1}(\mathbf{f}^{n+1} - \mathbf{K}\mathbf{w}^{n+1}) \\
\dot{\mathbf{w}}_{cd}^{n+1} &= \dot{\mathbf{w}}^n + \frac{1}{2}\Delta t (\ddot{\mathbf{w}}^{n+1} + \dot{\mathbf{w}}^n)
\end{aligned}
\quad (10) $$

For convenience, we summarize the front-shock inducing method presented in the preceding section as follows:

$$ 
\begin{aligned}
\mathbf{w}^{n+c} &= \mathbf{w}^n + \Delta t_c \dot{\mathbf{w}}^n + \frac{\Delta t_c^2}{2} \ddot{\mathbf{w}}^n \\
\ddot{\mathbf{w}}^{n+c} &= \mathbf{M}^{-1}(\mathbf{f}^{n+c} - \mathbf{K}\mathbf{w}^{n+c}) \\
\mathbf{w}_{fs}^{n+1} &= \mathbf{w}^n + \Delta t_c \alpha \dot{\mathbf{w}}^n + \Delta t_c^2 \beta_1(\alpha) \ddot{\mathbf{w}}^n + \Delta t_c^2 \beta_2(\alpha) \ddot{\mathbf{w}}^{n+c} \\
\ddot{\mathbf{w}}_{fs}^{n+1} &= \mathbf{M}^{-1}(\mathbf{f}^{n+1} - \mathbf{K}\mathbf{w}^{n+1}) \\
\dot{\mathbf{w}}_{fs}^{n+1} &= \dot{\mathbf{w}}^n + \frac{1}{2}\Delta t \{\ddot{\mathbf{w}}^n + \dot{\mathbf{w}}^{n+1}\} \\
\beta_1(\alpha) &= \frac{1}{6}\alpha(1 + 3\alpha - \alpha^2), \quad \beta_2(\alpha) = \frac{1}{6}\alpha(\alpha^2 - 1), \quad \alpha = \frac{\Delta t}{\Delta t_c}
\end{aligned}
\quad (11) $$

Let us now combine the two methods (10) and (11) to obtain

$$ 
\begin{aligned}
\mathbf{w}^{n+1} &= \theta \mathbf{w}_{fs}^{n+1} + (1 - \theta) \mathbf{w}_{cd}^{n+1}, \quad 0 \le \theta \le 1 \\
\dot{\mathbf{w}}^{n+1} &= \theta \dot{\mathbf{w}}_{fs}^{n+1} + (1 - \theta) \dot{\mathbf{w}}_{cd}^{n+1} \\
\ddot{\mathbf{w}}^{n+1} &= \theta \ddot{\mathbf{w}}_{fs}^{n+1} + (1 - \theta) \ddot{\mathbf{w}}_{cd}^{n+1}
\end{aligned}
\quad (12) $$

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
<figure><img src="image_9.png" /><figcaption>Figure 8. Responses to box wave by the new method (12) with $\Delta t = 0.5\Delta x/c$ and $\theta = 0.5$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Responses to box wave by the central difference method (10) that triggers post-shock oscillations and the front-shock oscillation method (11) when integrated with $\Delta t = 0.5\Delta x/c$.</figcaption></figure>

herein, we have performed additional numerical experiments for all ranges of the Courant numbers, namely $0 < c\Delta t/\Delta x \le 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.

## 4. 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><img src="image_11.png" /><figcaption>Figure 10. Bar consisting of two heterogeneous materials.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Extrapolation of the displacements ($w_1(x, t)$, $w_2(x, t)$) at steps ($(n+c1), (n+c2)$) followed by interpolations at the next step $(n+1)$ for heterogeneous problems.</figcaption></figure>

First, we obtain the pushforward displacements for the two domains, $(\mathbf{w}_1^{n+c1}, \mathbf{w}_2^{n+c1})$ as follows (See Fig. 11)

$$ 
\begin{aligned}
\mathbf{w}_1^{n+c1} &= \mathbf{w}_1^n + \Delta t_{c1} \dot{\mathbf{w}}_1^n + \frac{\Delta t_{c1}^2}{2} \ddot{\mathbf{w}}_1^n \\
\mathbf{w}_2^{n+c1} &= \mathbf{w}_2^n + \Delta t_{c1} \dot{\mathbf{w}}_2^n + \frac{\Delta t_{c1}^2}{2} \ddot{\mathbf{w}}_2^n
\end{aligned}
\quad (13) $$

It should be noted that $\Delta t_{c1}$ is the maximum pushforward timestep without triggering instability. Using these pushforward displacements, we obtain the accelerations $(\ddot{\mathbf{w}}_1^{n+c1}, \ddot{\mathbf{w}}_2^{n+c1})$ from the following partitioned equation set (see, e.g., [42] for details in the partitioned solution procedures):

$$ 
\begin{aligned}
\ddot{\mathbf{w}}_1^{n+c1} &= \mathbf{M}_1^{-1}(\mathbf{f}_1^{n+c1} - \mathbf{K}_1 \mathbf{w}_1^{n+c1} - \mathbf{B}_1 \boldsymbol{\lambda}_1^{n+c1}) \\
\ddot{\mathbf{w}}_2^{n+c1} &= \mathbf{M}_2^{-1}(\mathbf{f}_2^{n+c1} - \mathbf{K}_2 \mathbf{w}_2^{n+c1} - \mathbf{B}_2 \boldsymbol{\lambda}_2^{n+c1}) \\
\ddot{\mathbf{w}}_{2b}^{n+c1} - \mathbf{L}_{f1} \ddot{\mathbf{w}}_f^{n+c1} &= \mathbf{0}, \quad \mathbf{w}_{1b} = \mathbf{B}_1^T \ddot{\mathbf{w}}_2^{n+c1} \\
\ddot{\mathbf{w}}_{2b}^{n+c1} - \mathbf{L}_{f2} \ddot{\mathbf{w}}_f^{n+c1} &= \mathbf{0}, \quad \mathbf{w}_{2b} = \mathbf{B}_1^T \ddot{\mathbf{w}}_1^{n+c1} \\
\mathbf{L}_{f1}^T \boldsymbol{\lambda}_1^{n+c1} + \mathbf{L}_{f2}^T \boldsymbol{\lambda}_2^{n+c1} &= \mathbf{0}
\end{aligned}
\quad (14) $$

where $(\mathbf{B}_1, \mathbf{B}_2)$ are the Boolean matrices that extract the interface accelerations, $(\mathbf{L}_{f1}$ and $\mathbf{L}_{f2})$ are the extractor for the frame acceleration $(\ddot{\mathbf{w}}_f)$; $(\boldsymbol{\lambda}_1$ and $\boldsymbol{\lambda}_2)$ are the interface forces, respectively, as shown in Figure (12); and $(\mathbf{w}_{1b}$ and $\mathbf{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 $(\mathbf{w}_{1b}, \mathbf{w}_{2b})$ are related to the frame displacement $(\mathbf{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 $(\boldsymbol{\lambda}_1, \boldsymbol{\lambda}_2, \ddot{\mathbf{w}}_f)$ require, for the one-dimensional bar, the solution of a $(3 \times 3)$ matrix equation, and for multi-dimensional problems, a $(3n_g \times 3n_g)$ matrix

<figure><img src="image_13.png" /><figcaption>Figure 12. Heterogeneous bar partitioned. Here, we treat the displacement at the interface node, which is preserved as a frame displacement $(\mathbf{w}_f)$ and two distinct Lagrange multipliers are introduced to enforce the two independent interface compatibility conditions (cf., Equation 14).</figcaption></figure>
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 + c2)$ is obtained by

$$ 
\begin{aligned}
\mathbf{w}_2^{n+c2} &= \mathbf{w}_2^n + \Delta t_{c2} \dot{\mathbf{w}}_2^n + \frac{\Delta t_{c2}^2}{2} \ddot{\mathbf{w}}_2^n \\
\ddot{\mathbf{w}}_2^{n+c2} &= \mathbf{M}_2^{-1}(\mathbf{f}_2^{n+c2} - \mathbf{K}_2 \mathbf{w}_2^{n+c2} - \mathbf{B}_2 \lambda_2^{n+c2}), \quad \lambda_2^{n+c2} \approx \lambda_2^{n+c1}
\end{aligned}
\quad (15) $$

The pullback displacements at the time step $(n + 1)$ are thus obtained by

$$ 
\begin{aligned}
\mathbf{w}_{fs1}^{n+1} &= \mathbf{w}_1^n + \Delta t_{c1} \alpha \dot{\mathbf{w}}_1^n + \Delta t_{c1}^2 \beta_1(\alpha_1) \ddot{\mathbf{w}}_1^n + \Delta t_{c1}^2 \beta_2(\alpha_1) \ddot{\mathbf{w}}_1^{n+c1} \\
\mathbf{w}_{fs2}^{n+1} &= \mathbf{w}_2^n + \Delta t_{c2} \alpha \dot{\mathbf{w}}_2^n + \Delta t_{c2}^2 \beta_1(\alpha_2) \ddot{\mathbf{w}}_2^n + \Delta t_{c2}^2 \beta_2(\alpha_2) \ddot{\mathbf{w}}_2^{n+c2} \\
\beta_1(\alpha_1) &= \frac{1}{6}\alpha_1(1 + 3\alpha_1 - \alpha_1^2), \quad \beta_2(\alpha_1) = \frac{1}{6}\alpha_1(\alpha_1^2 - 1), \quad \alpha_1 = \frac{\Delta t}{\Delta t_{c1}} \\
\beta_1(\alpha_2) &= \frac{1}{6}\alpha_2(1 + 3\alpha_2 - \alpha_2^2), \quad \beta_2(\alpha_2) = \frac{1}{6}\alpha_2(\alpha_2^2 - 1), \quad \alpha_2 = \frac{\Delta t}{\Delta t_{c2}}
\end{aligned}
\quad (16) $$

In these equations, $(\Delta t_{c1}, \Delta 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:

$$ 
\begin{aligned}
\dot{\mathbf{p}}_1^{n+1} &= (\mathbf{f}_1^{n+1} - \mathbf{K}_1 \mathbf{w}_{fs1}^{n+1}) \\
\dot{\mathbf{p}}_2^{n+1} &= (\mathbf{f}_2^{n+1} - \mathbf{K}_2 \mathbf{w}_{fs2}^{n+1}) \\
\dot{\mathbf{p}}_g^{n+1} &= \mathbf{L}_g^T \left\{ \begin{array}{c} \dot{\mathbf{p}}_1^{n+1} \\ \dot{\mathbf{p}}_2^{n+1} \end{array} \right\} \\
\ddot{\mathbf{w}}_{gfs}^{n+1} &= \mathbf{M}_g^{-1} \dot{\mathbf{p}}_g^{n+1} \\
\left\{ \begin{array}{c} \ddot{\mathbf{w}}_{fs1}^{n+1} \\ \ddot{\mathbf{w}}_{fs2}^{n+1} \end{array} \right\} &= \mathbf{L}_g \ddot{\mathbf{w}}_{gfs}^{n+1}
\end{aligned}
\quad (17) $$

where $\mathbf{p}_{1,2}$ are the momentum vector of two domains, $\mathbf{L}_g$ is a Boolean operator that assembles the two domain displacements into the corresponding assembled displacement, and $\ddot{\mathbf{w}}_g$ is the assembled acceleration vector.

Third, we update the velocity for each domain via

$$ 
\begin{aligned}
\dot{\mathbf{w}}_{fs1}^{n+1} &= \dot{\mathbf{w}}_1^n + \frac{1}{2} \Delta t \{ \ddot{\mathbf{w}}_1^n + \ddot{\mathbf{w}}_{fs1}^{n+1} \} \\
\dot{\mathbf{w}}_{fs2}^{n+1} &= \dot{\mathbf{w}}_2^n + \frac{1}{2} \Delta t \{ \ddot{\mathbf{w}}_2^n + \ddot{\mathbf{w}}_{fs2}^{n+1} \}
\end{aligned}
\quad (18) $$

Fourth, we obtain the assembled velocity and displacement by

$$ \dot{\mathbf{w}}_{gfs} = \mathbf{L}_g^+ \left\{ \begin{array}{c} \dot{\mathbf{w}}_{fs1}^{n+1} \\ \dot{\mathbf{w}}_{fs2}^{n+1} \end{array} \right\}, \quad \mathbf{w}_{gfs} = \mathbf{L}_g^+ \left\{ \begin{array}{c} \mathbf{w}_{fs1}^{n+1} \\ \mathbf{w}_{fs2}^{n+1} \end{array} \right\} \quad (19) $$

where $\mathbf{L}_g^+$ is a pseudo-inverse of $\mathbf{L}_g$.

Fifth, we compute the post-shock oscillating response via

$$ 
\begin{aligned}
\mathbf{w}_{gp}^{n+1} &= \mathbf{w}_g^n + \Delta t \dot{\mathbf{w}}_g^n + \frac{\Delta t^2}{2} \ddot{\mathbf{w}}_g^n \\
\ddot{\mathbf{w}}_{gp}^{n+1} &= \mathbf{M}_g^{-1}(\mathbf{f}_g^{n+1} - \mathbf{K}_g \mathbf{w}_g^{n+1}) \\
\dot{\mathbf{w}}_{gp}^{n+1} &= \dot{\mathbf{w}}_g^n + \frac{1}{2} \Delta t (\ddot{\mathbf{w}}_g^{n+1} + \ddot{\mathbf{w}}_g^n)
\end{aligned}
\quad (20) $$

where the subscript $g$ denotes the assembled quantities.
Sixth, use (12) to obtain the acceleration, velocity and the displacement, which we repeat for completeness:

$$
\begin{align}
\mathbf{w}_g^{n+1} &= \theta \mathbf{w}_{gs}^{n+1} + (1-\theta) \mathbf{w}_{gp}^{n+1}, \quad 0 \le \theta \le 1 \notag \\
\dot{\mathbf{w}}_g^{n+1} &= \theta \dot{\mathbf{w}}_{gs}^{n+1} + (1-\theta) \dot{\mathbf{w}}_{gp}^{n+1} \tag{21} \\
\ddot{\mathbf{w}}_g^{n+1} &= \theta \ddot{\mathbf{w}}_{gs}^{n+1} + (1-\theta) \ddot{\mathbf{w}}_{gp}^{n+1} \notag
\end{align}
$$

For partitioned computations of the two heterogeneous bars, we obtain the partitioned accelerations, velocity and displacement by

$$
\left\{ \begin{matrix} \ddot{\mathbf{w}}_1^{n+1} \\ \ddot{\mathbf{w}}_2^{n+1} \end{matrix} \right\} = \mathbf{L}_g \ddot{\mathbf{w}}_g^{n+1}, \quad \left\{ \begin{matrix} \dot{\mathbf{w}}_1^{n+1} \\ \dot{\mathbf{w}}_2^{n+1} \end{matrix} \right\} = \mathbf{L}_g \dot{\mathbf{w}}_g^{n+1}, \quad \left\{ \begin{matrix} \mathbf{w}_1^{n+1} \\ \mathbf{w}_2^{n+1} \end{matrix} \right\} = \mathbf{L}_g \mathbf{w}_g^{n+1} \quad (22)
$$

This completes the computational sequence for a time step.

5. STABILITY AND ACCURACY ANALYSIS OF THE PRESENT ALGORITHM

It turns out that, for linear and linearized equations of motion, the stability of the present algo-
rithm 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):

$$
\begin{equation}
\begin{split}
& (\mathbf{w}^{n+1} - 2\mathbf{w}^n + \mathbf{w}^{n-1}) \\
& \quad + \left[ \mathbf{I} - \frac{1}{2}\alpha^2 \Delta t_c^2 \mathbf{M}^{-1} \mathbf{K} + \frac{1}{2}\theta \beta_2(\alpha) \Delta t_c^4 (\mathbf{M}^{-1} \mathbf{K})^2 \right] (\mathbf{w}^n - \mathbf{w}^{n-1}) \\
& \quad - \Delta t_c^2 [\alpha - \theta \beta_2(\alpha) \Delta t_c^2 \mathbf{M}^{-1} \mathbf{K}] \frac{\alpha}{2} \mathbf{M}^{-1} \mathbf{K} (\mathbf{w}^n + \mathbf{w}^{n-1}) = 0
\end{split}
\tag{23}
\end{equation}
$$

Numerical stability of the mentioned difference equation can be assessed by focusing on the
highest mode of the equations of motion. This can be obtained by introducing

$$
\mathbf{w} = \mathbf{T}\mathbf{q}, \quad \mathbf{T}^T\mathbf{M}\mathbf{T} = \mathbf{I}, \quad \mathbf{T}^T\mathbf{K}\mathbf{T} = \mathbf{\Lambda}, \quad \operatorname{diag}(\mathbf{\Lambda}) = \{\omega_1^2, \omega_2^2, \dots, \omega_N^2\} \quad (24)
$$

where $\omega_k$ is the undamped frequency for the $k$th mode, $\mathbf{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 $\mathbf{w}^{n+1}$ and $\mathbf{w}^n$ in terms of $\mathbf{w}^{n-1}$ by

$$
\mathbf{w}^{n+1} = \mu^2 \mathbf{w}^{n-1}, \quad \mathbf{w}^n = \mu^2 \mathbf{w}^{n-1} \tag{25}
$$

yielding for the highest frequency the following characteristic equation:

$$
\begin{align}
& (\mu - 1)^2 + \frac{1}{2}[\alpha^2 \Omega_{\max}^2 - \theta \beta_2(\alpha) \Omega_{\max}^4](\mu - 1) \nonumber \\
& \quad + [\alpha^2 - \alpha \theta \beta_2(\alpha) \Omega_{\max}^2] \frac{\Omega_{\max}^2}{2}(\mu + 1) = 0, \quad \Omega_{\max}^2 = \Delta t_c^2 \omega_{\max}^2 \tag{26} \\
& \beta_2(\alpha) = \frac{\alpha}{6}(\alpha^2 - 1), \quad \alpha = \frac{\Delta t}{\Delta t_c} \nonumber
\end{align}
$$

For stability, the amplification factor $\mu$ must be bounded by

$$
|\mu| \le 1 \tag{27}
$$

which can be mapped onto the entire left-hand plane of the z-plane given by

$$
\mu = \frac{1 + z}{1 - z} \tag{28}
$$
which, when substituted into (26) leads to

$$ 
\begin{aligned}
& az^2 + bz + c = 0 \\
& a = 4 - \Omega^2[\alpha^2 - \theta\beta_2]\Omega^2 \\
& b = \Omega^2[\theta\beta_2(1 - \alpha)\Omega^2] \\
& c = \alpha(\alpha - \theta\beta_2\Omega^2)\Omega^2
\end{aligned} 
\qquad (29) $$

For stability: $a \ge 0$, $b \ge 0$, $c \ge 0$ via Routh-Hurwitz criterion

For $\alpha \le 1$, we observe with $\gamma = \frac{1}{2}$ and $\theta = \frac{1}{2}$ that the conditions of $b \ge 0$ and $c \ge 0$ are satisfied. This leaves to examine $a$:

$$ a = 4 - \alpha^2 \Omega_{\max}^2 + \theta \beta_2 \Omega_{\max}^4 \qquad (30) $$

Figure13 shows the stability of limits of the present method plotted for stable zones versus the stepsize ratio ($\alpha = \Delta t / \Delta t_c$). In that figure, when $\alpha = 1$ (viz., $\Delta t = \Delta t_c$), the stability limit is given by

$$ \Delta t_c \le 2/\omega_{\max} \qquad (31) $$

When the stepsize ($\Delta$) is reduced from the critical stepsize ($\Delta t_c$), the sufficient stable stepsize is given by

$$ \Delta t \le 2\alpha/\omega_{\max} \qquad (32) $$

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 ($\alpha = \Delta t / \Delta t_c = 0.5$). Whereas the central difference method ($\alpha = 1$) introduces no damping ($|\mu| = 1.0$), the present method does introduce when $\alpha < 1.0$. Nevertheless, the amount of numerical damping is less than 7% when $\omega_{\max} \Delta t = 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 $\alpha = 0.5$ versus $\omega_{\max} \Delta t$ as shown in Figure 15. Note that the frequency is increased for the central difference method ($\alpha = 1$) whereas the present method with $\alpha = 0.5$ decreases the frequency. The absolute error of the present

<figure><img src="image_14.png" /><figcaption>Figure 13. Stability of the present method for varying stepsize ratios ($\alpha = \Delta t / \Delta t_c$).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14. Amplification (or damping) factor of the present method versus stepsize.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15. Frequency error of the present method versus stepsize.</figcaption></figure>

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.

# 6. 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.

## 6.1. 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
<figure><img src="image_17.png" /><figcaption>Figure 16. Responses to box wave by the central difference and the variational explicit methods for heterogeneous bar. Box wave propagation along bar (left end is free, right end fixed, time = 85).</figcaption></figure>

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 = 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.

<figure><img src="image_18.png" /><figcaption>Figure 17. Responses to box wave by the new method for heterogeneous bar. Box wave propagation along bar (left end is free, right end fixed, time = 85).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18. Responses to box wave by the new method for heterogeneous bar. Box wave propagation along bar (left end is free, right end fixed, time = 340).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19. Responses to box wave by the new method for heterogeneous bar. Box wave propagation along bar (left end free, right end fixed, time = 1020).</figcaption></figure>

## 6.2. 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 $Cr = c_{stiff} \Delta t / \Delta x = 0.85$ is chosen, the corresponding Courant number for the soft bar becomes $Cr = 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 = 1105$ and $t = 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].
<figure><img src="image_21.png" /><figcaption>Figure 20. Responses to box wave by the new method for heterogeneous bar. Box wave propagation along bar (left end free, right end fixed, time = 425).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21. Responses to box wave by the new method for heterogeneous bar. Box wave propagation along bar (left end free right end fixed, (left) time = 1105 and (right) time = 1700).</figcaption></figure>

### 6.3. Box wave applied at soft-stiff-stiffer and stiffer-stiff-soft heterogeneous bars

We have applied to a heterogeneous 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

$$f = 1.0, \quad \text{if } 0 \le T \le (L/c)_{\text{left}}, \quad T = N\Delta t \quad (33)$$

where $(L_{\text{left}}, c_{\text{left}})$ are the bar length and the speed of propagation of the left end bar; and, $(T, \Delta t, N)$ are the time duration of applied constant force, $\Delta t$ is the integration step size, and $N$ is the number of integration steps, respectively.
<figure><img src="image_23.png" /><figcaption>Figure 22. Bar consisting of three heterogeneous materials.</figcaption></figure>

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 $Cr < 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 $Cr_{stiffer} = 0.8$, which makes that of the soft segment to be $Cr_{soft} = 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.

<figure><img src="image_24.png" /><figcaption>Figure 23. Response to box wave for a bar consisting of three heterogeneous materials ($dt = 0.80$, $t = 1100$).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure 24. Response to box wave for a bar consisting of three heterogeneous materials ($dt = 0.80, t = 1250$).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 25. Response to box wave for a bar consisting of three heterogeneous materials ($dt = 0.80, t = 200$).</figcaption></figure>

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.

## 7. 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 $c\Delta t/\Delta x < 1$. The present method consists of a combination of two explicit methods: the classical central difference
<figure><img src="image_27.png" /><figcaption>Figure 26. Response to box wave for a bar consisting of three heterogeneous materials ($dt = 0.80, t = 800$).</figcaption></figure>

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 $Cr = c\Delta t/\Delta 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.

## ACKNOWLEDGEMENTS

The first author (KCP) has been partially supported by World Class University (WCU) Program through the Korea Science and Engineering Foundation funded by the Ministry of Education, Science and Technology, Korea (Grant Number R31-2008-000-10045-0). The second and third authors have been partially funded by the Brain Korea 21 Program (KAIST Valufacture Institute of Mechanical Engineering), Ministry of Education, Science and Technology.
REFERENCES

1. Courant R, Friedrichs K, Lewy H. On the partial difference equations of mathematical physics. *IBM Journal, March 1967*, pp. 215-234. English translation of the 1928 German original, "R. Courant, K. Friedrichs and H. Lewy, *ber die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen* 1928; **100**(1):32-74.

2. Lax P. Hyperbolic systems of conservation laws, II. *Communications on Pure and Applied Mathematics* 1957; **10**:537.

3. Lax PD, Wendroff B. Systems of conservation laws. *Communications on Pure and Applied Mathematics* 1960; **13**(2):217-237.

4. Richtmyer RD, Morton KW. *Difference Methods for Initial-Value Problems*, 2nd edn. Wiley-Interscience: New York, 1967.

5. Cherukuri HP. Dispersion analysis of numerical approximations to plane wave motions in an isotropic elastic solid. *Computational Mechanics* 2000; **25**(4):317-328.

6. Krenk S. Dispersion-corrected explicit integration of the wave equation. *Computer Methods in Applied Mechanics and Engineering* 2001; **191**:975-987.

7. Liu Y, Sen MK. A new timespace domain high-order finite-difference method for the acoustic wave equation. *Journal of Computational Physics* 2009; **228**:8779-8806.

8. Marfurt KJ. Accuracy of finite difference and finite element modeling of the scalar and elastic wave equation. *Geophysics* 1984; **49**:533-549.

9. Zahradnik J, Moczo P, Hron T. Testing four elastic finite-difference schemes for behavior at discontinuities. *Bulletin of the Seismological Society of America* 1993; **83**:107-129.

10. Zuazua E. Propagation, observation, and control of waves approximated by finite difference methods. *SIAM REVIEW* 2005; **47**(2):197-243.

11. Zuazua E. Optimal and approximate control of finite difference approximation schemes for the 1D wave equation. *Rendiconti di Matematica e delle sue Applicazioni* 2004; **24**:201-237.

12. Glimm J, Li XL, Liu Y, Xu Z, Zhao N. Conservative front tracking with improved accuracy. *SIAM Journal on Numerical Analysis* 2003; **41**(5):1926-1947.

13. Gloth O, Hanel D, Tran L, Vilsmeier R. A front tracking method on unstructured grids. *Computers & Fluids* 2003; **32**:547-570.

14. Leveque RJ, Shyue KM. Two-dimensional front tracking based on high resolution wave propagation methods. *Journal of Computational Physics* 1996; **123**:354-361.

15. Mao D-K. Toward front tracking based on conservation in two-dimensional space. *SIAM Journal on Scientific Computing* 2000; **22**(1):113-151.

16. Nguyen V-T, Peraire J, Khoo BC, Persson P-O. A discontinuous Galerkin front tracking method for two-phase flows with surface tension. *Computers & Fluids* 2010; **39**:1-14.

17. Becache E, Joly P, Tsogka C. An analysis of new mixed finite elements for the approximation of wave propagation problems. *SIAM Journal on Numerical Analysis* 2000; **37**:1175-1201.

18. Codina R. Finite element approximation of the hyperbolic wave equation in mixed form. *Computer Methods in Applied Mechanics and Engineering* 2008; **197**:1305-1322.

19. Guddati MN, Yue B. Modified integration rules for reducing dispersion error in finite element method. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**:275-287.

20. Harari I, Haham S. Improved finite element methods for elastic waves. *Computer Methods in Applied Mechanics and Engineering* 1998; **166**:143-164.

21. Harari I. Special issue on new computational methods for wave propagation. *Wave Motion* 2004; **39**(4):279-280.

22. Park KC. Practical aspects of numerical time integration. *Computers and Structures* 1977; **7**:343-353.

23. Park KC, Underwood PG. A variable-step central difference method for structural dynamic analysis, I: theoretical aspects. *Computer Methods in Applied Mechanics and Engineering* 1980; **22**:241-258.

24. Seriani G, Oliveira SP. Dispersion analysis of spectral element methods for elastic wave propagation. *Wave Motion* 2008; **45**(6):729-744.

25. TAMMA KK, NAMBURU RR. A new finite element based Lax-Wendroff/Taylor-Galerkin methodology for computational dynamics. *Computer Methods in Applied Mechanics and Engineering* 1988; **71**:137-150.

26. Yue B, Guddati MN. Dispersion-reducing finite elements for transient acoustics. *Journal of the Acoustical Society of America* 2005; **118**(4):2132-2141.

27. Huang H, Costanzo F. On the use of space-time finite elements in the solution of elasto-dynamic problems with strain discontinuities. *Computer Methods in Applied Mechanics and Engineering* 2002; **191**:5315-5343.

28. Hulbert GM, Hughes TRJ. Space-time finite element methods for second-order hyperbolic equations. *Computer Methods in Applied Mechanics and Engineering* 1990; **84**:327-347.

29. Mullen R, Belytschko T. Dispersion analysis of finite element semi-discretizations of the two-dimensional wave equation. *International Journal for Numerical Methods in Engineering* 1982; **18**:11-29.

30. Ainsworth M, Monk P, Muniz W. Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation. *Journal of Scientific Computing* 2006; **27**(1-3).

31. Li X, Yao D, Lewis RW. A discontinuous Galerkin finite element method for dynamic and wave propagation problems in non-linear solids and saturated porous media. *International Journal of Numerical Methods in Engineering* 2003; **57**:1775-1800.
32. Mabssout M, Pastor M. A Taylor-Galerkin algorithm for shock wave propagation and strain localization failure of viscoplastic continua. *Computer Methods in Applied Mechanics and Engineering* 2003; **192**:955-971.

33. Xing Y, Shu C-W. High-order well-balanced finite volume WENO schemes and discontinuous Galerkin methods for a class of hyperbolic systems with source terms. *Journal of Computational Physics* 2006; **214**:567-598.

34. Idesman A. Accurate time integration of linear elastodynamics problems. *CMES* 2011; **71**(2):111-148.

35. Smith JM. Recent developments in numerical integration. *Journal of Dynamic Systems, Measurement and Control* 96, Ser. G-1 1974; **1**:61-70.

36. Cho SS, Huh H, Park KC. A time-discontinuous implicit variational integrator for stress wave propagation analysis in solids. *Computer Methods in Applied Mechanics and Engineering* 2011; **200**:649-664.

37. Dehghana M, Saadatmandi A. Variational iteration method for solving the wave equation subject to an integral conservation condition. *Chaos, Solitons and Fractals* 2009; **41**:1448-1453.

38. Kane C, Marsden JE, Ortiz M, West M. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. *International Journal of Numerical Methods in Engineering* 2000; **49**:1295-1325.

39. Lew A, Marsden JE, Ortiz M, West M. Variational time integrators. *International Journal of Numerical Methods in Engineering* 2004; **60**:153-212.

40. Wendlandt JM, Marsden JE. Mechanical integrators derived from a discrete variational principle. *Physica D* 1997; **106**:223-246.

41. Newmark NM. A method of computation for structural dynamics. *Journal of Engineering Mechanics, Proceedings of ASCE* 1959; **85**:67-94.

42. Park KC, Felippa CA. A variational principle for the formulation of partitioned structural systems. *International Journal of Numerical Methods in Engineering* 2000; **47**:395-418.

43. Dharmaraja S, Kesari H, Darve E, Lew A. Time integrators based on approximate discontinuous Hamiltonians. *Proceedings of the 3rd International Conference on Computer Methods in Structural Dynamics & Earthquake Engineer*, Corfu, Greece, 2011; 25-28.

44. Romero I, Orden JCG. Applications of thermodynamically consistent time integrators in solids, structural, and multi-body dynamics. *Proceedings of the 3rd International Conference on Computer Methods in Structural Dynamics & Earthquake Engineer*, Corfu, Greece, 2011; 29-31.