High-Order Eulerian Simulations of Multimaterial Elastic–Plastic Flow

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 [7–9], 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 [12–14] 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.

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.

where υ₀ = 10⁻³² is used to prevent division by zero in smooth regions.

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.

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

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

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.

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⁻⁶, 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.

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⁻⁶ 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.

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 γ₁ = 2, p_∞,1 = 68.23 GPa, μ₁ = 39.38 GPa, ρ_0,1 = 8930 kg/m³, and σ_Y,1 = 0.12 GPa. The EOS parameters for aluminum are γ₂ = 2.088, p_∞,2 = 32.98 GPa, μ₂ = 27.09 GPa, ρ_0,2 = 2712 kg/m³, 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 γ₁p_∞,1 and density scale ρ_0,1.

where p₁ = 5 × 10⁻² is the nondimensional preshock pressure and p₂ 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 T_osc is quantified and plotted as a function of κ^−1∕2. A linear trend between T_osc and κ^−1∕2 is observed in Fig. 4 as predicted by the theory.

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 η₀k = 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.

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

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 Y_Cu > 0.5, and the location of bubbles is defined as the minimum x for which Y_Cu < 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.

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.

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 η₀k. 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 η₀k. For η₀k of 0.1, the mixing width remains small. For larger η₀k, the deposited vorticity renders the interface unstable, and leads to rapid growth of the mixing width especially after reshock. Increasing η₀k 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.

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.

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.

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.

• •

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

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.