An Efficient B-Spline Lagrangian/Eulerian Method for Compressible Flow, Shock Waves, and Fracturing Solids

An extended amount of mechanical literature exists on the traditional partitioned methods for coupling compressible flow with solids. The most well-known categories include the arbitrary Lagrangian Eulerian method [Banks et al. 2016], the immersed boundary method [Pasquariello et al. 2016; Wang et al. 2017; Monasse et al. 2012], the discontinuous Galerkin method [Kosík et al. 2015], the ghost fluid method [Bailoor et al. 2017], and the Lattice Boltzmann method [Li and Favier 2017]. Many such methods only support small deformation or vibration. Some of them support large deformation without fracture. In computer graphics, Sewall et al. [2009] also adopted the partitioned method and weakly coupled shockwaves with rigid bodies.

A characteristic visual feature of explosions and shockwaves is their destructive effects on structures. Without modeling fracture mechanics, the above works cannot simulate destructive visual effects. Furthermore, these methods’ partitioned (instead of monolithic) approaches often suffer from numerical instabilities when large deformation occurs.

To enable fracturing of solids, we adopt the recently advanced Material Point Method (MPM) [Sulsky et al. 1995; Stomakhin et al. 2013; Zhang et al. 2016; Jiang et al. 2016] for the solid domain. Nowadays, MPM is widely used in computer graphics and computational engineering for simulating elastoplastic materials with topology change. Due to the hybrid usage of meshless Lagrangian particles and Eulerian grids, MPM automatically captures material separation and contact.

MPM also naturally simulates mixed materials [Jiang et al. 2016; Hu et al. 2018] without the need for modeling the explicit coupling among them. However, due to the usage of a single grid, MPM suffers from artificial adhesion artifacts at material interfaces. Furthermore, it is computationally expensive and memory intensive to represent gas with particles. Therefore, we only model the solid phase with MPM and adopt an Eulerian grid-based discretization for the fluid phase. Existing work has coupled MPM sediments with Eulerian incompressible fluids [Gao et al. 2018a; Baumgarten et al. 2021]. Ma et al. [2009]’s method switches from MPM to FDM during the detonation process based on whether the material interaction or the gas pressure evolution is dominant. However, they do not model the coupling between the two phases.

In summary, most existing schemes for coupling compressible flow with deformable solids do not simultaneously support strong coupling, solid fracture, and efficient integration of the compressible flow. The motivation of our work is thus to develop a new computationally efficient, monolithic, and fracture-supporting strong coupling scheme for destructive effects due to shockwaves and high-speed flows.

A highly related work is the Interface Quadrature Material Point Method (IQ-MPM) by Fang et al. [2020]. IQ-MPM uses MPM particles to represent both solid and incompressible liquid phases. Strong coupling is achieved by modeling the interfacial interaction with a layer of imaginary quadrature points. Their approach is computationally efficient, only requiring the solution of an SPD system per time step. We extend some of their ideas from incompressible liquids to compressible flow.

The last piece of our framework concerns the performance of the compressible flow solver. We follow Kwatra et al. [2009];, 2010]’s work, which successfully performed the time-splitting method in the compressible flow regime, removed the strict restriction of the sound speed and built an SPD pressure-only system—a system that can be integrated into the coupling phase from [Fang et al. 2020]. In recent years, the semi-implicit solver in [Kwatra et al. 2009;, 2010] was further improved in its stability [Grétarsson et al. 2011], mass conservation [Grétarsson and Fedkiw 2013], and subgrid resolution [Hyde and Fedkiw 2019]. These work couple gas with sediments and simple elastic bodies. They also capture simple fracture events by plugging pressure-induced forces into a rigid body fracture tool. On the other hand, our work supports arbitrarily complex fracture events of nonlinear elastic solids under extreme deformation, as well as splitting and merging dynamics of plastic materials, without extra efforts beyond having an MPM solver for modeling these processes.

In summary, we present a new method for simulating potentially destructive interactions between compressible flow (shock waves) and nonlinear elastoplastic solids. Our technical contributions include

In contrast to Kwatra et al. [2010]’s method, which couples staggered Marker-and-Cell (MAC) grid Finite Difference-based fluids with purely Lagrangian solids, our method is based on the non-staggered grid, which avoids extrapolation in the advection step for the compressible flow. Our framework is based on operator splitting, and it remains stable under large time step sizes determined by the pressure gradient and the maximum velocity. All of our examples ran with a CFL number 0.5.

Following standard literature, here, we describe our governing equations: a conservative form of the inviscid compressible Euler equations for fluids and the conservation of mass and momentum for elastoplastic solids. Introducing density \(\rho\) , velocity \(\mathbf {u}\) , and pressure p, the compressible flow ( \(\bullet ^f\) ) region is governed byand the solid region ( \(\bullet ^s\) ) is governed byAt the solid-gas interface we enforce the slip velocity condition and pressure continuity:In Equation (3), E is the fluid’s enthalpy (total energy per unit volume). It relates to the internal energy e throughThe system is closed by the equation of state (EOS) for an ideal gas [Kwatra et al. 2010]where we pick \(\gamma =1.4\) .

In the remaining subsections, we split these coupled equations into advection and non-advection parts. For the reasons mentioned earlier in the introduction, we discretize the gas with an Eulerian grid and the solid with the hybrid Lagrangian/Eulerian MPM.

We summarize the overall pipeline in Figure 7. Our method consists of the following procedures in each time step.

Firing bullets. We use high-temperature high-pressure gas to fire bullets at varying conditions in Figure 5. A heavy bullet reaches the transonic state. It does not break the wavefront until near the end of the simulation and creates uniform vortices behind its tail. The lighter bullet reaches the supersonic state. It breaks the wavefront early in the simulation and produces massive chaotic vortices. In another test, we add a silencer at the end of the firing chamber to observe its successful suppression of the primary shock. This bullet has its speed slightly decreased, and it moves forward, passing the much slower shockwave and slicing out a secondary shock. Pumpkin explosion. We simulate the explosions of a plastic pumpkin and an elastic pumpkin by setting up ellipsoidal high-pressure gas inside them (Figure 6). They are both severely torn apart. The deformation of the plastic pumpkin is permanent, while the elastic pumpkin jiggles vividly in the back-flow. Supernova. We use the planar wave to mimic a far-field supernova explosion, destroying a mini solar system in Figure 8. Due to the plastic nature of the solids in this example, they are easily blown into dust. Cylinder flow and von Kármán vortices. We simulate the interactions between March 3 flow and cylinders made of different elastoplastic materials. Varying the Young’s modulus and switching plasticity on and off, we observe intricate fluid-solid coupling behaviours and distinct patterns of von Kármán vortices; see Figure 10. Airplane in a wind tunnel. We put a toy airplane in a wind tunnel and gradually increase the flow speed until fractures occur; see Figure 11. The loss of balance causes the plane to spin around in the tunnel. The plane is severely deformed and fractured during its dynamic interaction with the supersonic flow. Cactus near explosion. We simulate the interaction between a cactus forest and a supersonic wind blow caused by a far-field explosion in Figure 12. The cactus plants are cut off near their roots as the shock approaches. Mach Diamond and Koalas. The Mach diamond appears at the exhaust of a slightly over-expanded jet with an outlet pressure a bit lower than the ambient. The jet is compressed and expands periodically to form the shape of diamonds (shown in Figure 18). In a 3D simulation (Figure 19), we put elastic koala toys with different densities above air jets. We pick the densities of the koalas so that the lightest koala is easily blown up, the heaviest one breaks the wavefront, and the middle one stays balanced for a while before slipping off the exhaust area.

The solid advection takes place during the grid-to-particle transfer at the end of the previous time step [Jiang et al. 2016]. At the beginning of the current step, the solid particles transfer their quantities to the Eulerian grid via a standard particle-to-grid (P2G) operation. Subsequently, the location and the velocity of the interface between the two phases are obtained.

The fluid advection refers to the process of updating the conserved fluid variables. They are treated as independent variables, temporally ignoring the influence of the pressure. The fluid is updated to an intermediate state ( \(\bullet ^*\) ) after the advection. Here, we use the finite difference method and the second-order WENO-LLF flux reconstruction to solve the advection equations:

In simulations involving massive destruction, many solid fragments are scattered all over. The resulting moving interface configurations between the two phases are complex, and require specific attention to avoid errors. We discuss our approaches below.

Mixed reflective and passable interface. The P2G operation may result in cells with small solid volume fractions. The advection at the interface between such a cell and a gas cell is neither reflective nor passable. Hence, we take the following three steps to blend between purely reflective and passable interfaces. Firstly, we model every interface as a mixture of a reflective wall and a passable passage (see Figure 3 in [Hyde and Fedkiw 2019]). The passable ratio R is derived from the cross-section of the solid’s volume by treating it as a ball. Thus, we havewhere \(V^s\) is the solid’s volume [Baumgarten and Kamrin 2019; Gao et al. 2018a]. Secondly, we extrapolate both the conserved variables and the velocity for the passable flux stencil, and extrapolate only the conserved variables, but mirror reflect the velocity for the reflective flux stencil [Forrer and Jeltsch 1998; Forrer and Berger 1999]. Lastly, the obtained two fluxes are blended using the passable ratio. We illustrate this mixed interface \(I_1\) in Figure 9.

Avoiding using reflection at any interface but the one being evaluated flux at. Consider an interface such as \(I_1\) in Figure 9, there are scenarios where the adjacent gas cell \(C_1\) is touching another solid cell \(C_0\) . In such cases, the ghost values of the higher-order terms in the flux stencil (those in \(C_0\) ) become ambiguous. In particular, it becomes a question of whether to treat the interface \(I_0\) as a reflective interface. If we choose to do so, this ambiguity becomes more severe in multi-dimensional cases since there may exist other solid-gas interfaces in other dimensions (for example, at \(C_0\) ’s top face). In such scenarios, reflecting the ghost value at different interfaces yields inconsistent results and causes simulation artifacts. To avoid this issue, we do not perform reflection at any solid-gas interfaces other than the one we are evaluating the flux at. We use the extrapolated value from the nearest gas cell corresponding to a lower-order term in the stencil for the additional required higher-order terms from solid cells.

After the advection step, we update the fluid’s pressure with the EOS using the intermediate quantities (instead of using semi-Lagrangian methods as in [Kwatra et al. 2009]) to avoid the numerical drifting caused by non-conservative advection methods [Aanjaneya et al. 2013].

Projection is a common term in the simulation of incompressible flow [Bridson 2015] and refers to projecting the intermediate post-advection velocities onto a divergence-free subspace, using pressure to enforce the divergence-free constraints. In the compressible flow regime, velocity and pressure variables do not appear in the advection step as we use the conserved variables. Nonetheless, we treat velocity and pressure variables as primary states and describe the equations they must obey in the following.

In summary, through operator splitting, we can conduct a post-advection coupled solve for solid-fluid coupling. The coupling step occurs after the solid’s particle-to-grid transfer and the fluid’s advection, but before the MPM grid update (see Figure 7). More specifically, the time-discretized equations are

Combined with the interfacial continuity Equations (6) and (7) at the solid-fluid interface \(\Gamma\) and the Dirichlet wall boundary conditions for \(\mathbf {u}^{\bullet }\) at \(\partial \Omega ^{\bullet }\) , these equations form a first-order system in terms of the unknown fields \(\lbrace \mathbf {u}^{f,n+1}, p^{f,n+1}, \mathbf {u}^{s,n+1}, p^{g,n+1}\rbrace\) and can be further spatially discretized with the finite element method.

Using B-spline shape functions for finite elements has been recently shown effective for incompressible flow [Fang et al. 2020; Gagniere et al. 2020]. Here, we extend the formulation to our compressible fluid-structure coupled system.

For continuity of elastic forces in MPM, the solid’s velocity field must be \(\mathcal {C}^1\) continuous [Hughes 2012; Bardenhagen and Kober 2004] (requiring the shape function to be quadratic polynomials), while the pressure can be either piecewise linear or piecewise constant. Fang et al. [2020] described two choices for discretizing coupled incompressible fluids and solids. The first scheme adopts quadratic B-splines for velocities and linear functions for pressures on both phases, resulting in a B2B1-B1-B2B1 scheme, where the B1 in the middle refers to the discretization order for the interfacial pressure. To reduce memory bandwidth and computational cost, they also described a second scheme B2B0-B0-B1B0, switching to piecewise constant pressures everywhere and piecewise linear velocities for the fluid phase.

In our case, however, if B2B0-B0-B1B0 is applied, the velocity nodes will be staggered for the two phases, and the interface will cut through the inner volume of the fluid cell. This “cut-cell” will require extra tedious treatments to the advection step in the fluid’s boundary cell due to our usage of the flux-based solver.

To mitigate this problem, we stick to B2B1 for the solids and add offsets to the fluid coordinates, taking advantage of the Eulerian nature of the fluid. In most of our examples, the solid phase only occupies a small portion of the whole domain; thus, the use of B2B1 does not severely affect the performance. We refer to our discretization choice as B2B1-B1/B0-B1B0 and use it in large-scale 3D simulations. We use B2B1 for the fluid phase in our 1D and 2D benchmarks.

The weak form of our coupled system in Section 2.5 can then be obtained by properly applying test functions \(\mathbf {q}^{f,s}, r^{f,s}\) :

The weak form of the velocity conditions at the solid-fluid interface and the domain boundaries areHere, \(\mathbf {q}_{\alpha }\) collocates with the velocity nodes and is nonzero only on its \(_{\alpha }\) component, \(\Gamma\) is the solid-fluid interface, p is the pressure, \(\mathbf {u}_b\) is the normal velocity at the slip boundaries of the solid or fluid domains, and r is the test function collocated with the pressure nodes. For B2B1, we would havewhere \(N^k\) is the kth order B-spline kernel and \(\delta\) is the Kronecker delta. A discrete and compact form of the above equations can be obtained through the standard FEM assembling process, which yieldswhere the \(\mathbf {h}\) and \(\mathbf {Y^{f,s}}\) are the stacked vectors of the pressure DOFs at the solid-fluid interface and the domain boundaries, respectively, which appear because of the integration by parts of \(\int _{\Omega } \nabla p^{n+1} \cdot \mathbf {q}_{\alpha } d\Omega\) . The right hand side is

Equation (34) as a KKT system is challenging to solve. We transform it into an SPD system by mass lumping and eliminating the velocity DOFs. The mass matrix \(\mathbf {M}^{f,s}\) and the stiffness matrix \(\mathbf {S}^{f,s}\) are not ill-conditioned. Nevertheless, the lumping makes them diagonal-scaling matrices and easily invertible. We eliminate velocity DOFs using

This yields an SPD and well-conditioned pressure-only system:Here,and \(rhs=(\mathbf {S}^{f,-1} \mathbf {P}^{f,*}+ \Delta t \mathbf {G}^{f,T} \mathbf {U}^{f,*}_{\alpha }, \Delta t \mathbf {B}^{f,T} \mathbf {U}^{f,*}_{\alpha } - \Delta t b^{f}, \Delta t \mathbf {H}^{f,T} \mathbf {U}^{f,*}_{\alpha } + \Delta t \mathbf {H}^{s,T} \mathbf {U}^{s,*}_{\alpha }, \Delta t \mathbf {B}^{s,T} \mathbf {U}^{s,*}_{\alpha } - \Delta t b^{s}, \mathbf {S}^{s,-1} \mathbf {P}^{s,*} + \Delta t \mathbf {G}^{s,T} \mathbf {U}^{s,*}_{\alpha })^T\) . This system is diagonal dominant and can be solved efficiently with Jacobi-preconditioned Conjugate Gradient (CG). It usually takes less than 80 CG iterations to converge for our examples.

After solving the linear system, we use Equation (36) to update the solid’s velocities and use local gradients to update the conserved variables of the fluid [Kwatra et al. 2009]. We found in our experiments that updating the solid’s velocities with local gradients sometimes causes instability, possibly owing to the nature of the quadratic order of the velocity discretization (see Figure 13 for a visual comparison). On the other hand, we cannot use Equation (35) to update the fluid’s velocities because otherwise we cannot update the total energy properly.

The velocities on the solid cells after the coupling step have taken the contributions from the fluid and require a final standard MPM grid update step to incorporate hyperelasticity. Note that Fang et al. [2020] adopted a different order and put the MPM step before the coupling step. Doing that, however, results in instabilities in our framework when a strong shock interacts with stiff materials (see Figure 14 for a visual comparison).

Traditional one-step approaches for simulating compressible flow use the characteristic decomposition to handle inlet and outlet boundary conditions. In time-splitting approaches, Kwatra et al. [2009] elaborated on the stencil modification for calculating boundary fluxes in the advection step. We need to properly model the Neumann boundary condition in our coupled system specifying the pressure gradient in Equation (14). Equivalently we specify the final \(t^{n+1}\) velocities at the boundary faces and modify the right hand side of Equation (37). The user specifies the inlet’s velocity values directly, and the outlet takes the intermediate velocity values after advection.

In summary, to model inlets and outlets, we

1D validation. We validate the convergence behavior of our method following the 1D cases from [Kwatra et al. 2009]. We generate high resolution (4,096 cells) results using their solver and measure the squared density deviation asThe results are plotted in Figure 16, showing that our method has excellent convergence behaviors for relatively higher speed flows, while at lower speeds the convergence order is 1.

2D validation. We use the central circular shock test from [Kwatra et al. 2009] as a 2D validation. The results at varying resolutions (ranging from 256 to 2048) are compared against a 4,096-resolution benchmark obtained by their method. We also use this test to perform an ablation study comparing the accuracy between B2B1 and B1B0 velocity/pressure orders. The results (see Figure 16) indicate that B2B1 and B1B0 both converge under refinement, showing unnoticeable differences in the visual results. We thus adopt B1B0 in our large-scale examples for its lower computational cost.

Boundary treatment. We validate our treatment for inlet and outlet boundaries using the classic stair flow test from computational fluid dynamics. The result in Figure 15 is visually plausible.

Fluid-solid interaction. We validate the qualitative accuracy of our fluid-solid coupling scheme using the classic lifting cylinder test in 2D. The cylinder is sampled with uniformly distributed material points with a high Young’s modulus. Our result is shown in Figure 17. The cylinder is lifted to the top of the channel toward the outlet, reproducing the result in [Forrer and Berger 1999].

We show more examples “(refer to Table 1 for detailed parameters)” to demonstrate the efficacy of our method with an emphasis on the natural advantage of MPM’s handling of large deformation and topology change.

Car driving through a land mine. We simulate a car driving through a triggered landmine in Figure 1. The car’s main body is plastic, while its wheels are made of much lighter elastic materials. The intense explosion tears apart and blows away the rear wheels and causes permanent distortion to the body.

“Tribunny”. We simulate three bunny toys surrounding an explosion in Figure 2. The heavy elastic bunny and the plastic bunny are damaged without being pushed far away. On the other hand, the lighter elastic bunny is fully blown away by the shockwave, as expected. Elastic shell explosion. We simulate an exploding elastic shell in Figure 3. The shell is filled with highly compressed hot gas. The gas expands and breaks the shell. The process first results in two layers of primary shockwaves spreading outwards. Later on, it triggers an inward injection of the gas and the mixing of the post-explosion exhaust. The elastic fragments naturally interact with the complex flow pattern. Firing Cannon at a truck. We simulate projectiles hitting truck toys in Figure 4. The trucks are first damaged by the projectiles and then deformed more by the shockwaves. The material parameters of each truck determine the severity of its final damage.

Our method contains several limitations. Firstly, even though our mixed grid treatment increases subgrid resolution to some extent, it cannot account for detailed sub-grid solid geometry. For example, our method cannot capture the accurate reflection at a sub-grid non-axis-aligned solid surface. Cut cell methods may improve the accuracy here but will also bring in more challenges.

Secondly, the computational cost of our method imposes restrictions on the practical resolution we can achieve. The performance bottleneck of our method is the construction of the linear system due to its requirement of performing many sparse matrix multiplications. In our implementation, we accelerate this part based on pre-obtained knowledge of the sparsity pattern. Our optimization strategy cannot be easily generalized to more sophisticated future extensions using cut-cells or adaptive meshes.

A third limitation comes from the shortcomings of operator splitting. In fluid-structure interactions, the velocity at the interface greatly impacts the flow pattern. Like almost every other method using operator splitting, our method does not maintain the same interface velocity at different stages of the algorithm. Our work does not cover the impact on the accuracy of such inconsistency.

Lastly, we simplify the treatment of the inlet/outlet boundary conditions by directly specifying the desired velocities. A more accurate treatment of them in a time-splitting framework necessitates more rigorous mathematical descriptions of the coupling phase from the characteristic decomposition of the flux terms.

In this work, we rely on MPM’s artificial numerical fracture for simulating the fragmentation of solids. It is an exciting future work to incorporate continuum damage mechanics [Wolper et al. 2019;, 2020] into the MPM component of our framework to model more versatile brittle and ductile fracture mechanics.

It would also be interesting to couple shockwaves with incompressible liquids, granular media, and soil sediments. For powder-like solid particles, our mixed reflective and passable interface scheme in the advection step may no longer fully capture the turbulent sub-grid details; thus, new schemes are required.

From the performance perspective, we observe that while the Adaptive Mesh Refinement (AMR) schemes are well developed in the traditional characteristic-based compressible flow solvers, they are not broadly studied in time-splitting methods. It is thus meaningful to design a spatially adaptive treatment for the coupled phases, extending Gao et al. [2017].

As mentioned in limitations, the performance of our solver is severely bottlenecked by the linear system’s construction. It is a promising future work to devise matrix-free solvers with multigrid preconditioners [Wang et al. 2020a] and mitigate the solver to modern computational platforms following [Gao et al. 2018; Wang et al. 2020b].