We develop a new high-order numerical method for continuum simulations of multimaterial phenomena in solids exhibiting elastic–plastic behavior using the diffuse interface numerical approximation. This numerical method extends an earlier single material high-order formulation that uses a tenth-order high-resolution compact finite difference scheme in conjunction with a localized artificial diffusivity (LAD) method for shock and contact discontinuity capturing. The LAD method is extended here to the multimaterial formulation and is shown to perform well for problems involving shock waves, material interfaces and interactions between the two. Accuracy of the proposed approach in terms of formal order (eighth-order) and numerical resolution is demonstrated using a suite of test problems containing smooth solutions. Finally, the Richtmyer–Meshkov (RM) instability between copper and aluminum is simulated in two-dimensional (2D) and a parametric study is performed to assess the effect of initial perturbation amplitude and yield stress.

## Introduction

Solids undergoing large deformations in contact with fluids occur in several engineering contexts including several manufacturing processes and processes involving cavitation. Robust and accurate handling of the interface between materials is the primary challenge encountered in numerical simulations of such phenomena. Lagrangian or arbitrary Lagrangian–Eulerian methods, where the mesh distorts with the motion of the interface, often fail when dealing with problems involving large deformations. Eulerian methods, on the other hand, are attractive since they employ fixed grids that do not deform with the motion and are not susceptible to failure due to mesh entanglement. They are also ideally suited for problems involving wave propagation and fluid-like behavior. However, tracking the interface between materials becomes more complicated and requires a more careful treatment.

The motion of a single solid undergoing nonlinear, elastic–plastic deformations can be tracked using the fully Eulerian framework introduced by Godunov and Romenskii [1,2] and by Plohr and Sharp [3,4]. Several improvements to the formulation were suggested in Refs. [5] and [6]. These latter formulations rely on the inverse deformation gradient tensor as the fundamental kinematic variable for tracking deformations of the solid. The constitutive behavior of solids is described by hyperelastic equations of state that maintain thermodynamic consistency. Extensions of this Eulerian framework to multiple materials follow either the sharp-interface methodology [79], or the diffuse-interface approximation [10,11]. The sharp-interface approach is most often used with level-set methods [8,9]. Such an approach keeps the interface sharp and Riemann solvers or ghost cells for the interface jump conditions can be used with relative ease. However, these methods leak mass since they are nonconservative and also tend to be excessively dissipative at the interface. Conservative level set methods that reduce the magnitude of conservation errors have been proposed [1214] but they are not discretely conservative and have not been used previously in the context of elastic–plastic behavior in solids. Diffuse-interface approaches, on the other hand, do not keep the interface sharp, but numerically smear it over a few grid points. They also require an explicit mixture treatment since there are a few numerically mixed cells close to the interface. These are, however, conservative and can be used with low dissipation schemes unlike the sharp-interface approach.

In addition to handling material interfaces, the capability to handle shocks and contact discontinuities is required in simulations of compressible solids and fluids. Most previous work has focused on implicitly adding the required dissipation in order to capture shocks and material interfaces. Ndanou et al. [11] employ a Godunov-type scheme for numerical solution of the diffuse-interface multimaterial formulation, which is first-order accurate. The hybrid central difference—weighted essentially non-oscillatory scheme employed by López Ortega et al. [9], while formally third-order accurate, introduces upwinding and significant numerical dissipation. An alternative methodology for capturing shocks and interfaces in fluids (mainly ideal gases) was proposed and developed by Cook [15,16] and Kawai et al. [17], and combines high-resolution central compact finite difference schemes [18] with explicitly added numerical dissipation. In our recent work [19], this localized artificial diffusivity (LAD) scheme was extended and shown to be successful in capturing shocks in a single compressible material spanning a large range of continuum models (from gases, liquids, to elastic–plastic solids) in a unified manner. In the present work, high-order compact finite differences in conjunction with the LAD scheme are further extended to problems involving interactions of shocks with material interfaces, particularly between solids undergoing elastic–plastic deformations. This high-order method is inspired by the use of similar high-order compact finite differences in multimaterial problems in gases (see Refs. [20] and [21]), where the benefits of high resolution of compact finite differences are significant compared to lower order schemes. Such high-order compact finite difference schemes were shown to be competitive with lower-order schemes such as fifth-order weighted essentially non-oscillatory schemes from the point of view of computational cost per grid point as well (see Table 1 in Ref. [21]).

The mathematical formulation and numerical solution procedure are discussed in Sec. 2. The diffuse-interface multimaterial model is used, since it assumes that a mixture of materials exists at every point in the computational domain, and is suitable for use with global numerical schemes such as high-order compact finite differences. One-dimensional (1D) test problems that establish order of accuracy of the code are discussed, followed by a two-dimensional (2D) demonstration problem involving interaction between a shock and a perturbed interface between copper and aluminum. The paper ends with a brief summary and comments about future work.

## Governing Equations and Numerical Methodology

In this section, we briefly describe the governing equations and numerical formulation used in this work focusing mainly on the nuances associated with the multimaterial aspects. For a detailed discussion on the single material formulation and numerical method, we refer the reader to Ref. [19]. We use the convention that bold symbols denote vector variables and underlined symbols denote second-order tensor variables.

The fundamental equations governing motion of a multimaterial continuum in Eulerian form are the mass, momentum, and energy equations shown in index notation below:
$∂ρYm∂t+∂ρYmuk∂xk=−∂(Jm∗)i∂xi$
(1)
$∂ρui∂t+∂∂xk[ρuiuk−σik]=∂τik∗∂xk$
(2)
$∂∂t[ρ(ε+12ujuj)]+∂∂xk[ρuk(ε+12ujuj)−σikui]=∂∂xk[τik∗ui−qk∗]$
(3)
where $u$ is the velocity, and Ym is the mass fraction of material m. ρ, ε, and $σ¯$ are the mixture density, internal energy, and Cauchy stress, respectively, given in terms of species-specific quantities ρm, εm and $σ¯m$, by
$ρ=∑m=1Mαmρm, ε=∑m=1MYmεm, σ¯=∑m=1Mαmσ¯m$
(4)
where αm is the volume fraction of material m and M is the total number of materials. The volume fractions have to satisfy the constraint $∑m=1Mαm=1$. $τ¯∗=2μ∗S¯+(β∗−(2/3)μ∗)(∇·u)I¯$ is the artificial viscous stress and $q∗=−κ∗∇T+∑m=1MhmJm∗$ is the artificial heat flux. $S¯=(1/2)(∇u+(∇u)T)$ is the strain rate tensor, T is the temperature, and μ*, β*, and κ* are the artificial shear viscosity, artificial bulk viscosity, and artificial thermal conductivity, respectively. The form of μ*, β,* and κ* used in this paper is the same as in Ref. [19]. The second term in the expression for $q∗$ is the enthalpy diffusion term [16] where hm = εm + pm/ρm is the enthalpy of species m. $(Jm∗)i=−ρ(Dm∗(∂Ym/∂xi)−Ym∑kDk∗(∂Yk/∂xi))$ is the Fickian diffusive flux for species m in direction xi and $Dm∗$ is the artificial diffusivity of species m. The form for $Dm∗$ is similar to Ref. [16] and is given by
$Dm∗=max{CDcs|∑k=13∂rYm∂xkrΔkrΔYm2|¯,CYcs2(|Ym|−1+|1−Ym|)ΔYm¯}$
(5)
where cs is the speed of sound (or linear longitudinal wave speed) in the mixture, Δk is the grid spacing in the xk direction, r = 4 is a parameter that localizes the dissipation to just the interfacial region, and the overbar $(f¯)$ denotes a truncated Gaussian filter applied along each grid direction to smooth out sharp cusps introduced by the absolute value operator. Details on the truncated Gaussian filter are given in Ref. [16]. CD and CY are model constants that control the interface thickness and the mass fraction over- and under-shoot, respectively. CD = 3 × 10−3 and CY = 100 were found to yield robust results for all the problems considered here while keeping the mass fraction over- and under-shoot to within 2–3% and the numerical interface thickness to ∼ 4–5 grid points. $ΔYm$ is a weighted grid spacing parameter given by
$ΔYm=∑i=13Δi|∂Ym∂xi|∑k=13(∂Ym∂xk)2+υ0$
(6)

where υ0 = 10−32 is used to prevent division by zero in smooth regions.

To track the deformations in the solid medium, we solve for only the elastic part of the inverse deformation gradient tensor $g¯e$, and model the effect of plasticity. The equations governing the evolution of the elastic part of the inverse deformation gradient tensor for material m, $g¯me$ can be written as in Refs. [9] and [11]
$∂(gme)ij∂t+∂uk(gme)ik∂xj=uk(∂(gme)ik∂xj−∂(gme)ij∂xk)+ζg(ρmρ0,m|g¯m|−1)(gme)ij+12μm(ρmρ0,m)τ0,m[R(||σ¯m′||2−23σY,m2)μm2](gme)ikσ′kj$
(7)

The second term on the right-hand side is to ensure consistency with the material density and ζg = 1/(6Δt) is a time-step dependent constant based on stability considerations. The last term on the right-hand side of Eq. (7) accounts for plasticity effects where μm denotes the shear modulus, ρ0,m, density in the undeformed state, and $σ¯m′$ is the deviatoric part of $σ¯m, σY,m$ is the yield stress, and 1/τ0,m is an inverse relaxation timescale of material m. The Ramp function, R(), turns off plasticity effects whenever the yield criterion $||σ¯m′||2=σ¯m′:σ¯m′=σij,m′σij,m′<23σY,m2$.

In this work, only perfectly plastic solid behavior is considered. Furthermore, in the test problems considered here, we assume τ0,m/τf ≪ 1, where τf is a flow time scale measure so that the materials relax to the yield surface almost instantaneously. This assumption is justified for the materials considered in this work (copper and aluminum) but for materials where the plastic relaxation timescales are comparable to the flow timescales, the plastic source terms can be treated explicitly.

### Pressure and Temperature Equilibrium.

In addition to the Eqs. (1)(3) and (7), we require two more equations per material apart from the equation of state (EOS) to close the system. Equations for material compaction (described by the volume fraction) and species energy can be derived based on thermodynamic considerations [22] to close the system. Kapila et al. [23] reduce the full governing equations to reduced ones based on asymptotic analysis and small timescale arguments for pressure and temperature relaxation processes. In this work, we assume pressure and temperature equilibrium between different materials present within one grid cell. This assumption is justified if the pressure and temperature relaxation timescales are much smaller than the timescales resolved in the continuum model. This is typically true for materials with similar properties like two metals or two gases but is not true in general for vastly dissimilar materials as demonstrated in Ref. [23]. With this assumption, for a total of M species, we have
$p=p1=p2=⋯=pM$
(8)
$T=T1=T2=⋯=TM$
(9)
$∑m=1Mαm=1, ∑m=1MYmεm=ε, εm=εm(ρm,pm,g¯me)$
(10)
which is a system of 2M + 2 nonlinear equations that is solved for the 2M + 2 unknowns, p, T, αm and εm after each stage of the time-advancement scheme. This method is an extension of a similar procedure described in Ref. [16] for ideal gases, to the current system involving elastic–plastic solids.

#### Constitutive Description.

An EOS is required to close the above system of equations. In this study, we consider a classical hyperelastic form of the EOS, without additional contributions due to capillary/Korteweg effects. The stress is derived from an energy functional based on thermodynamic consistency through the Clausius–Duhem inequality. The internal energy and the resulting Cauchy stress of material m are given by
$εm=pm+γmp∞,m(γm−1)ρm+μm4ρ0,mtr((g¯me−1¯)2)$
(11)
$σ¯m=−pm1¯−μmρmρ0,m(|G¯m|−2/3(G¯m2)′−|G¯m|−1/3(G¯m)′)$
(12)

where $G¯m=(g¯me)Tg¯me$ is the elastic Finger tensor, $(A¯)′$ denotes the deviatoric part of the tensor $A¯$, and $1¯$ is the identity tensor. Material behavior is described by EOS parameters, γm, p,m, shear modulus, μm, and the undeformed density, ρ0,m. This constitutive description is the same as Refs. [11] and [19].

### Numerical Solution Procedure.

We solve the set of Eqs. (1)(7) on a fixed Eulerian structured grid using a tenth-order compact finite difference scheme [18]. A fourth-order five stage Runge–Kutta time integrator [24] is used. The stiff plastic source terms are integrated implicitly, following the method detailed in Ref. [25]. After each stage, the conserved variables are filtered for de-aliasing using an eighth-order compact filter [18]. Compact finite difference schemes provide spectral-like resolution capabilities but are much more flexible than spectral schemes. The centered schemes used here have no inherent dissipation. Dissipation required for shock and interface capturing is added explicitly through the artificial viscous stress, artificial conductive flux, and artificial diffusive flux in Eqs. (1)(3) as described in Refs. [16] and [19].

For the correct limiting behavior when the volume fractions and mass fractions simultaneously approach zero, the species density of species m is regularized using the density estimate obtained from the inverse deformation gradient tensor and is given by
$ρm=ρYm+υ0ρ0,m|g¯me|αm+υ0$
(13)

## Results

Results from verification studies for the single-material formulation that shows the high-order and high-resolution ability of the proposed numerical method are shown in Ref. [19]. In this section, we focus on results for the multimaterial formulation starting with some one-dimensional verification tests followed by the problem of shock interaction with a copper–aluminum material interface.

### Order of Accuracy.

In this section, we present two test problems to quantify the accuracy of the method in capturing material interfaces. For problems involving shock waves, shear discontinuities and entropy waves in a single material, we refer the reader to Ref. [19].

In order to establish the order of accuracy of the proposed numerical method, we consider advection of a smooth interface between two elastic solids on a periodic domain. The problem is nondimensionalized based on a length scale L*, a pressure scale, p* and a density scale, ρ* that gives a velocity scale $u∗=p∗/ρ∗$. We take the domain length equal to L* so that the domain in nondimensional units is x ∈ [0, 1]. The two solids have identical properties in nondimensional units $γ1=γ2=1.4, p∞,1=p∞,2=1, μ1=μ2=1, σY,1=σY,2=∞$, except for the undeformed densities, which are ρ0,1 = 1, and ρ0,2 = 100. These quantities have been nondimensionalized with scales $p∗=p∞,1∗=p∞,2∗$ and $ρ∗=ρ0,1∗$ where $p∞,m∗$ and $ρ0,m∗$ are the dimensional values of EOS parameter p and the undeformed density of material m. The volume fraction of material 1 varies from αmin to 1 − αmin, αmin = 10−6, with a Gaussian profile initially centered at x = 0.5 and with standard deviation 0.1. Both materials are initialized with dimensionless velocity u = 0.5. Due to the advection of the material mixture the Gaussian pulse is seen to be centered at x = 0.75 at t = 0.5 in Figs. 1(a) and 1(b). Using the analytical solution of the advecting mixture, a convergence study is performed. Figure 1(c) shows that the errors in volume and mass fractions and mixture density reduce at eighth-order.

Fig. 1
Fig. 1
Close modal

The second example considered consists of a linear plane strain wave of sufficiently small strength crossing an interface between two materials with density ratio of 2. The material constants are $γ1=γ2=1.4, p∞,1=p∞,2=1, μ1=μ2=0, ρ0,1=0.5$, and ρ0,2 = 1. The reference scales used here for nondimensionalization are $p∗=p∞,1∗=p∞,2∗$, and $ρ∗=ρ0,2∗$. A slab of material 2 between x = 0.5 and 0.9 is surrounded by material 1, and the domain in x ∈ [0, 1] is periodic. Each material interface is described by volume fractions varying from αmin = 10−6 to 1 − αmin with error function profiles of thickness 0.001. A compression wave is initialized at x = 0.35 with normal stress profile $σ11=−10−4 exp[−((x−cL,1t−0.35)/0.035)2]$ where $cL,i=(λL,i+2μL,i)/ρ0,i$ is the linear longitudinal wave speed in material i and μL,i = μi and λL,i = γip,i − 2μi/3 are the equivalent Lamé coefficients (see, e.g., Ref. [26]). Other components of stress, deformations, densities, and velocities are initialized so as to obtain a linear plane strain wave propagating to the right.

Interaction of the Gaussian pulse with the interface at x = 0.5 leads to a transmitted as well as a reflected pulse, as seen in Figs. 2(a) and 2(b). The transmission and reflection coefficients, which determine the magnitudes of the respective waves, are given analytically in the case of a sharp interface by Achenbach [26]. Error convergence rates are analyzed in Fig. 2(c). In a first set of calculations, the interface thickness is set to be 1/NX so that the diffuse–interface solution approaches the sharp-interface limit with increasing number of grid points. Figure 2(c) indicates that reducing the interface thickness with the grid resolution (so that number of grid points in the interface is kept fixed across grid resolutions) leads to convergence to the sharp-interface analytical solution at second-order. Since the numerically diffuse interface converges to the sharp interface at only first-order, degradation of the order of convergence is expected here. Nevertheless, it is encouraging that the convergence is still higher than first-order that would be expected. In a second set of calculations, the interface thickness is kept fixed equal to 0.01. Taking a very finely resolved (NX = 12,800) numerical solution as the reference solution, the errors reduce at eighth-order, consistent with the numerical schemes used.

Fig. 2
Fig. 2
Close modal

### Richtmyer–Meshkov Instability of a Copper–Aluminum Interface.

A two-dimensional test involving interaction of a shock with an interface between copper (Cu) and aluminum (Al) is discussed. Copper is described by the material parameters γ1 = 2, p,1 = 68.23 GPa, μ1 = 39.38 GPa, ρ0,1 = 8930 kg/m3, and σY,1 = 0.12 GPa. The EOS parameters for aluminum are γ2 = 2.088, p,2 = 32.98 GPa, μ2 = 27.09 GPa, ρ0,2 = 2712 kg/m3, and σY,2 = 0.297 GPa. These parameter values are derived from the Romenskii EOS parameters used for these metals and reported in Ref. [27]. The parameter values used here ensure that the speed of sound of a linear wave is equal to that with the Romenskii EOS used in Ref. [27]. All results are nondimensionalized by a pressure scale γ1p,1 and density scale ρ0,1.

#### Comparison With Linear Theory.

In this section, the Cu–Al Richtmyer–Meshkov problem described above is considered but with plasticity effects excluded in order to be able to compare to linearized theory of Plohr and Plohr [28]. The rectangular simulation domain extends from the point (x, y) = (−2, 0) to (6, 1) for lengths 8 and 1 in the x and y directions, respectively. A numerical sponge layer, which absorbs outgoing waves and minimizes reflection is applied at the left and right boundaries in the regions x ∈ [−2, 0] and x ∈ [4,6], respectively (see the Appendix for details). The top and bottom boundaries are periodic. All simulations are performed on a uniform and isotropic grid. The interface is initialized at x = 2 with a single-mode sinusoidal perturbation
$η(y)=2+η0 sin (ky)$
(14)
Since we use a diffuse interface approximation, we smear the initial interface over a few grid points. The initial volume fraction of Cu is given by
$αCu=12[1−(1−2αmin)erf(x−η(y)δtΔ)]$
(15)
where αmin = 10−6, δt = 3, Δ is the grid spacing, and erf() is the error function. Following Ref. [27], the parameter k = 4π is used so that two full wavelengths are included in the domain, and the initial amplitude of perturbations is given by η0k = 0.2. A shock moving to the right toward the Cu–Al interface with pressure ratio of 25 is initiated in copper at xs = 1.5 using the Rankine–Hugoniot relations derived from Eqs. (1)(3) and Eq. (7). Like the material interface, the shock is also numerically smeared over a few grid points. The pressure profile across the shock is given by
$p=p1+12(p2−p1)[1−erf(x−xsΔ)]$
(16)

where p1 = 5 × 10−2 is the nondimensional preshock pressure and p2 is the nondimensional postshock pressure. The velocity and density profiles across the shock are also similarly smoothed.

Plohr and Plohr [28] performed linearized analysis of the elastic RM problem and showed that for a heavy-light configuration as is the case here, the interface perturbation amplitude reduces after interaction with the shock and oscillates around a mean value with a characteristic frequency. As the shear moduli of the two materials is changed by a factor κ from μm to κμm, this characteristic frequency changes proportional to κ1∕2. The perturbation amplitude can be defined as half of the difference between the location of the spike and bubble. This method, however, is sensitive to the location of the interface on the Eulerian grid and the interpolation method used to find the location of the spike and bubble. Hence, we use an integral mixing width measure to find the perturbation amplitude $A(t)=4∫04Y¯CuY¯Aldx$ where $()¯$ indicates averaging in the transverse direction. This integral metric gives the perturbation amplitude for a linear or sawtooth interface exactly, but only up to a constant factor for sinusoidal interfaces. Figure 3 shows the interface perturbation amplitude normalized by its initial value as a function of time for different values of κ. We see that for all the curves, the interface amplitude reduces after shock interaction and oscillates with a characteristic frequency that increases with κ. In Fig. 4, the time period of this oscillation Tosc is quantified and plotted as a function of κ−1∕2. A linear trend between Tosc and κ−1∕2 is observed in Fig. 4 as predicted by the theory.

Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal

#### Plasticity Effects in the Copper–Aluminum RM Problem.

In this section, we consider the Cu–Al RM problem with plasticity effects. The simulation setup more closely resembles that in Ref. [27]. The rectangular simulation domain extends from the point (x, y) = (–2, 0) to (4, 1) for lengths 6 and 1 in the x and y directions, respectively. A portion of the domain is seen in Fig. 5. A numerical sponge layer, which absorbs outgoing waves and minimizes reflection, is applied at the left boundary in the region x ∈ [–2, 0], while a rigid wall is assumed at the right boundary at x = 4. The top and bottom boundaries are periodic. The interface is initialized at x = 2 with a single-mode sinusoidal perturbation with the same parameters as in Sec. 3.2.1, but with an initial amplitude given by η0k = 0.4. A shock with pressure ratio of 25 is initiated in copper at x = 1.5, and is moving to the right toward the Cu–Al interface. An x-t diagram of the problem with a flat interface is shown in Fig. 6(c). Different grid resolutions are labeled by the number of points in the transverse direction NY. A plot of the initial density profile (see Fig. 5) shows the location of the shock and perturbed material interface.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal

Results of a simulation using 768 × 128 grid points are seen in Figs. 6 and 7. The shock deposits vorticity on the material interface due to baroclinic torque. The deposition of baroclinic vorticity renders the interface unstable. The propagating shocks are strong enough that the yield stresses of both, aluminum and copper, are exceeded, and both solids undergo plastic flow. The vorticity contours at t = 1 (see Fig. 7(c)) show that, due to plastic flow, the bulk of the deposited vorticity is retained in the interfacial region, and only a small amount is carried away by weak elastic shear waves (not seen in the scale used for vorticity contours here). The density contours in Fig. 7(a) show that the interface is set in motion toward the right and has undergone a phase reversal since the shock interacts with a heavy-light interface. The shock transmitted into aluminum rebounds from the rigid wall at the right boundary at t ≈ 0.98. This reflected shock interacts with the disturbed interface a second time at t ≈ 1.4, leading to a left-moving transmitted shock and a right-moving, weaker, reflected shock. The reflected shock rebounds from the rigid wall once again, and a series of re-shock and interface interaction events is set up, eventually halting the translation of the interface. Since most of the vorticity is retained at the interface due to plastic effects and only a small portion is transported out by elastic shear waves, the interface curls up and a large deformation of the interface is observed as in Figs. 7(b) and 7(d) along with the characteristic shapes of the copper “spikes” and aluminum “bubbles.”

Fig. 7
Fig. 7
Close modal

The initiation of motion and linear and nonlinear growth of the instability are seen in the x locations of the spikes and bubbles in Fig. 6(a). The location of spikes is defined as the maximum x for which YCu > 0.5, and the location of bubbles is defined as the minimum x for which YCu < 0.5. The shock–interface interaction and ensuing interface instability is also seen in the evolution of integrated positive vorticity, $Γ+=∫(ωz+|ωz|/2)dxdy$, where $ωz=(∂u/∂y−∂v/∂x)$.

The effect of grid resolution on the density contours is shown in Fig. 8. As the grid is refined, more of the fine-scale structure of the interface is resolved. The location of shock waves in the domain is also seen to be the same across the different resolutions, which is a consequence of the fact that the formulation is fully conservative. Statistics plotted in Fig. 6 are well converged before the first reshock at t ≈ 1.4. After the first reshock, the differences between the NY = 64 and the NY = 128 simulation results are smaller than those between the NY = 32 and NY = 64 results. This indicates that simulations are approaching convergence, but are not fully converged on the 768 × 128 grid. We can quantify this further by using Richardson extrapolation to estimate the order of convergence. Just after reshock, at t = 0.375, we get an order of convergence of 2.82. Although not the eighth-order accuracy observed in smooth problems, the approximately third-order convergence is higher than the first-order convergence in approximating the sharp interface using a diffuse interface. However, after reshock at t = 1.4, we see that the order of accuracy is lower at 1.09, and the same first-order convergence in the circulation is observed after the second reshock as well. This degradation of the order of accuracy is because the fine-scale structure of the interface just before the first and second reshocks critically influences the vorticity deposited on the interface. Potentially, much higher grid resolutions are required in order to get into the asymptotic range of convergence post reshock, where convergence at a rate higher than first-order might be expected.

Fig. 8
Fig. 8
Close modal

The main benefit of using a high-order compact finite difference scheme is to get access to high numerical resolution over a large range of resolvable wavenumbers, similar to that afforded by a spectral method [18]. Although an exact, fair, comparison is not possible due to the different EOSs employed to describe copper and aluminum, the results in Ref. [27] (see Figure C.9) suggest that using their numerical method, effectively 1024 grid points in the y direction are required to resolve the features seen using the numerical method presented here with 64 points in the y direction (see Fig. 9(a)). This highlights the high resolution obtained with the present approach.

Fig. 9
Fig. 9
Close modal

#### Effect of Initial Amplitude.

The effect of changing the amplitude of the initial perturbation of the interface is seen in Fig. 9, which shows Cu mass fraction contours at time t = 4 for different values of η0k. Corresponding mixing width and net integrated positive vorticity are plotted in Fig. 10. The amount of vorticity deposited by the initial shock increases with increasing η0k. For η0k of 0.1, the mixing width remains small. For larger η0k, the deposited vorticity renders the interface unstable, and leads to rapid growth of the mixing width especially after reshock. Increasing η0k results in an increase in the baroclinic torque and hence the vorticity deposited on the interface as shown in Fig. 10(b). This leads to an increase in the extent of curl-up of the interface and inter-penetration of the spikes and bubbles.

Fig. 10
Fig. 10
Close modal

#### Effect of Yield Stress.

Figure 11 shows the effect of increasing the yield stress by a factor of κY compared to the baseline yield stress given in Sec. 3.2 for both materials. We see that increasing the yield stress inhibits the growth of the instability significantly. Figure 12 shows the density and vorticity contours at t = 0.5 and t = 1 in a purely elastic Cu–Al RM problem that corresponds to κY so that plastic effects are absent and only elastic effects are present. Comparing this to Figs. 7(a) and 7(c), we see that the main difference is that in the elastic case, the baroclinic vorticity deposited by the shock at the interface is transported away from the interface by elastic shear waves. For RM instability in fluids, baroclinic vorticity at the interface causes the highly nonlinear roll-up of the interface. Since vorticity is transported away from the interface in elastic materials, the interface is stabilized and no roll-up of the interface is observed. As plasticity effects are introduced by lowering the yield stress, vorticity transported away from the interface reduces since the magnitude of shear waves transporting vorticity out is lower. This subsequently enhances the growth of the interfacial instability.

Fig. 11
Fig. 11
Close modal
Fig. 12
Fig. 12
Close modal

## Conclusions

The main aim of this work is to enable high-order Eulerian simulations of large, elastic–plastic deformations of solids. To this end, we have developed a high-order numerical method composed of a tenth-order compact finite difference scheme and the LAD scheme for numerical regularization of shocks and interfaces to solve the governing equations for multimaterial elastic–plastic flows based on a diffuse interface approximation. The accuracy of the method was demonstrated through one-dimensional tests. The interaction of a shock wave with an interface separating copper and aluminum was simulated without plasticity effects and was shown to compare well with linear theory. The Cu–Al Richtmyer–Meshkov instability similar to that of Lopez-Ortega [27] was simulated and grid convergence results showed the superior resolution ability of the proposed numerical method. Further, the effect of initial amplitude was quantified. A study on the effect of yield stress showed that the RM instability is inhibited with increasing yield stress as a result of vorticity transported away from the interface by elastic shear waves. Future efforts will focus on Richtmyer–Meshkov flow problems in converging (cylindrical and spherical) geometries, and on problems involving solid–fluid interfaces with larger density ratios than those considered here. Modifications to the algorithm with the purpose of handling sliding at solid interfaces, and effects of thermal nonequilibrium, will also be explored in the future.

## Acknowledgment

We thank Dr. Andrew Cook for fruitful discussions, particularly related to the pressure and temperature equilibrium algorithm. We also thank the anonymous referees for their constructive comments.

## Funding Data

• Lawrence Livermore National Laboratory, Department of Energy (Grant No. B612155).

### Appendix: Details of the Numerical Sponge Layer

Details for the numerical sponge layer at the left x boundary used in Sec. 3.2 are presented here. The sponge region is defined by a sponge center location xs (−1.5 in Sec. 3.2) and a characteristic width δs. Dissipation is added in the sponge region x < xs by filtering the solution in the sponge region using a truncated Gaussian filter denoted by $()¯$ (see Ref. [16] for details on the Gaussian filter) and smoothly increasing the dissipation from zero in the interior. For a field f, the filtered field $f̃$ is given by
$f̃=(1−Θ)f+Θf¯$
(A1)

where $Θ=(1/2)[1−tanh(x−xs/δs)]$ is a blending function that is zero in the interior and one in the sponge region so that dissipation is only added in the sponge region. At each substep of the time integration scheme, this filtering procedure is applied to the primitive variables (pressure, velocities, density, volume fractions, and components of the inverse deformation gradient tensor) so as to dissipate any outgoing waves. In the problems described in this paper, the sponge was placed sufficiently far away so that the flow is essentially 1D. δs = 0.2 was found to be sufficient and absorbed any outgoing waves without reflections back into the domain.

## References

1.
Godunov
,
S. K.
, and
Romenskii
,
E. I.
,
1972
, “
Nonstationary Equations of Nonlinear Elasticity Theory in Eulerian Coordinates
,”
J. Appl. Mech. Tech. Phys.
,
13
(
6
), pp.
868
884
.
2.
Godunov
,
S. K.
, and
Romenskii
,
E. I.
,
2003
,
Elements of Continuum Mechanics and Conservation Laws
,
Springer
, New York.
3.
Plohr
,
B.
, and
Sharp
,
D.
,
1989
, “
A Conservative Eulerian Formulation of the Equations for Elastic Flow
,”
,
9
(4), pp.
481
499
.
4.
Plohr
,
B.
, and
Sharp
,
D.
,
1992
, “
A Conservative Formulation for Plasticity
,”
,
13
(4), pp.
462
493
.
5.
Trangenstein
,
J. A.
, and
Colella
,
P.
,
1991
, “
A Higher-Order Godunov Method for Modelling Finite Deformation in Elastic-Plastic Solids
,”
Commun. Pure Appl. Math.
,
44
(1), pp.
41
100
.
6.
Miller
,
G. H.
, and
Colella
,
P.
,
2001
, “
A High-Order Eulerian Godunov Method for Elastic-Plastic Flow in Solids
,”
J. Comput. Phys.
,
167
(1), pp.
131
176
.
7.
Hill
,
D. J.
,
Pullin
,
D. I.
,
Ortiz
,
M.
, and
Meiron
,
D. I.
,
2010
, “
An Eulerian Hybrid WENO Centered-Difference Solver for Elastic-Plastic Solids
,”
J. Comput. Phys.
,
229
(24), pp.
9053
9072
.
8.
Barton
,
P. T.
,
Deiterding
,
R.
,
Meiron
,
D. I.
, and
Pullin
,
D. I.
,
2013
, “
Eulerian Adaptive Finite-Difference Method for High-Velocity Impact and Penetration Problems
,”
J. Comput. Phys.
,
240
, pp.
76
99
.
9.
López Ortega
,
A.
,
Lombardini
,
M.
,
Pullin
,
D. I.
, and
Meiron
,
D. I.
,
2014
, “
Numerical Simulation of Elastic-Plastic Solid Mechanics Using an Eulerian Stretch Tensor Approach and HLLD Riemann Solver
,”
J. Comput. Phys.
,
257
(Part A), pp.
414
441
.
10.
Favrie
,
N.
, and
Gavrilyuk
,
S. L.
,
2012
, “
Diffuse Interface Model for Compressible Fluid—Compressible Elastic-Plastic Solid Interaction
,”
J. Comput. Phys.
,
231
(7), pp.
2695
2723
.
11.
Ndanou
,
S.
,
Favrie
,
N.
, and
Gavrilyuk
,
S. L.
,
2015
, “
Multi-Solid Multi-Fluid Diffuse Interface Model: Applications to Dynamic Fracture and Fragmentation
,”
J. Comput. Phys.
,
295
, pp.
523
555
.
12.
Olsson
,
E.
, and
Kreiss
,
G.
,
2005
, “
A Conservative Level Set Method for Two Phase Flow
,”
J. Comput. Phys.
,
210
(
1
), pp.
225
246
.
13.
Desjardins
,
O.
,
Moureau
,
V.
, and
Pitsch
,
H.
,
2008
, “
An Accurate Conservative Level Set/Ghost Fluid Method for Simulating Turbulent Atomization
,”
J. Comput. Phys.
,
227
(
18
), pp.
8395
8416
.
14.
Desjardins
,
O.
, and
Pitsch
,
H.
,
2009
, “
A Spectrally Refined Interface Approach for Simulating Multiphase Flows
,”
J. Comput. Phys.
,
228
(
5
), pp.
1658
1677
.
15.
Cook
,
A. W.
,
2007
, “
Artificial Fluid Properties for Large-Eddy Simulation of Compressible Turbulent Mixing
,”
Phys. Fluids
,
19
(5), p.
055103
.
16.
Cook
,
A. W.
,
2009
, “
Enthalpy Diffusion in Multicomponent Flows
,”
Phys. Fluids
,
21
, p.
055109
.
17.
Kawai
,
S.
,
Shankar
,
S. K.
, and
Lele
,
S. K.
,
2010
, “
Assessment of Localized Artificial Diffusivity Scheme for Large-Eddy Simulation of Compressible Turbulent Flows
,”
J. Comput. Phys.
,
229
(5), pp.
1739
1762
.
18.
Lele
,
S. K.
,
1992
, “
Compact Finite Difference Schemes With Spectral-Like Resolution
,”
J. Comput. Phys.
,
103
(1), pp.
16
42
.
19.
Ghaisas
,
N. S.
,
Subramaniam
,
A.
, and
Lele
,
S. K.
,
2016
, “
High-Order Eulerian Methods for Elastic-Plastic Flow in Solids and Coupling With Fluid Flows
,”
AIAA
Paper No. 2016-3350.
20.
Tritschler
,
V. K.
,
Olson
,
B. J.
,
Lele
,
S. K.
,
Hickel
,
S.
,
Hu
,
X. Y.
, and
,
N. A.
,
2014
, “
On the Richtmyer–Meshkov Instability Evolving From a Deterministic Multimode Planar Interface
,”
J. Fluid Mech.
,
755
, pp.
429
462
.
21.
Johnsen
,
E.
,
,
J.
,
Bhagatwala
,
A. V.
,
Cabot
,
W. H.
,
Moin
,
P.
,
Olson
,
B. J.
,
Rawat
,
P. S.
,
Shankar
,
S. K.
,
Sjögreen
,
B.
,
Yee
,
H.
, Zhong, X., and Lele, S. K.,
2010
, “
Assessment of High-Resolution Methods for Numerical Simulations of Compressible Turbulence With Shock Waves
,”
J. Comput. Phys.
,
229
(
4
), pp.
1213
1237
.
22.
Baer
,
M.
, and
Nunziato
,
J.
,
1986
, “
A Two-Phase Mixture Theory for the Deflagration-to-Detonation Transition (ddt) in Reactive Granular Materials
,”
Int. J. Multiphase Flow
,
12
(
6
), pp.
861
889
.
23.
Kapila
,
A.
,
Menikoff
,
R.
,
Bdzil
,
J. B.
,
Son
,
S.
, and
Stewart
,
D. S.
,
2001
, “
Two-Phase Modeling of Deflagration-to-Detonation Transition in Granular Materials: Reduced Equations
,”
Phys. Fluids
,
13
(
10
), pp.
3002
3024
.
24.
Kennedy
,
C. A.
,
Carpenter
,
M. H.
, and
Lewis
,
R. M.
,
2000
, “
Low-Storage, Explicit Runge–Kutta Schemes for the Compressible Navier–Stokes Equations
,”
Appl. Numer. Math.
,
35
(
3
), pp.
177
219
.
25.
Favrie
,
N.
, and
Gavrilyuk
,
S. L.
,
2010
, “
Dynamics of Shock Waves in Elastic-Plastic Solids
,”
ESAIM Proceedings
, Marseille, France, July 20–Aug. 28, pp.
1
18
.https://hal.inria.fr/file/index/docid/492411/filename/Nonlinear_plasticity_ESAIM-1.pdf
26.
Achenbach
,
J. D.
,
1975
,
Wave Propagation in Elastic Solids
,
Elsevier
, Amsterdam, The Netherlands.
27.
López Ortega
,
A.
,
2013
, “Simulation of Richtmyer-Meshkov Flows for Elastic-Plastic Solids in Planar and Converging Geometries Using an Eulerian Framework,”
Ph.D. thesis
, California Institute of Technology, Pasadena, CA.https://thesis.library.caltech.edu/7488/
28.
Plohr
,
J. N.
, and
Plohr
,
B. J.
,
2005
, “
Linearized Analysis of Richtmyer-Meshkov Flow for Elastic Materials
,”
J. Fluid Mech.
,
537
, pp.
55
89
.