Improved CE/SE scheme with particle level set method for numerical simulation of spall fracture due to high-velocity impact

Abstract

In the present paper, an Eulerian scheme combined with the hybrid particle level set method for numerical simulation of spall fracture due to high-velocity impact is proposed. An axisymmetric framework is established, based on an improved CE/SE scheme, to solve the high-velocity impact problems with large deformations, high strain rates and spall fractures. The hybrid particle level set method is adopted for tracking material interfaces and describing the formation and propagation of a crack. A novel representation of crack by level set is proposed. Numerical simulations are carried out and compared to the corresponding experimental results. The numerical results are in good agreement with the experimental data. The edge effects are reproduced and the decrease of scab thickness with increase in impact velocity is observed owing to the numerical analysis. It is proved that our computational technique is feasible and reliable for analyzing the spall fracture.

Introduction

Spall fracture happens mainly in the field of high-velocity impact. It is a kind of dynamic fracture that occurs inside a material when two decompression waves interact and produce enough tension to cause damage. The process of the spall fracture is considered as a sequence of nucleation, growth and coalescence of microscopic voids or cracks [1]. The classical plate-impact experiment is usually used to study the mechanisms of the spall fracture in the laboratories [2], [3], [4].

Numerical analysis plays an important role in the study of the spall fracture. It enables the researchers to obtain the information such as pressure, strain rate and porosity evolution inside the target during the tests, which cannot be captured during the experiments. Two kinds of dynamic fracture models for the numerical simulation of the spall, which consider the microvoiding process, are currently popular. They are the void growth (VG) model [5], [6], [7], [8], [9] and the nucleation and growth (NAG) model [3], [10]. The VG model is a mathematical model of ductile holes growth under the application of a mean tensile stress. The NAG model considers more physical details including the nucleation and the growth processes. Eftis et al. [11] performed numerical simulations, by their own constitutive-microdamage equations, for describing the impact of spherical projectiles with rectangular target plates at a speed of 6.0 km/s. Three ratios of the projectile diameter to the target thickness were chosen, and different damage features were presented in the simulations. Ikkurthi and Chaturvedi [12] adopted four damage models for the simulations of the spall fracture in metal plates. The best match with the experiments was obtained by the VG model. Several physical phenomena such as the edge effects were also studied. Czarnota et al. [10] carried out the finite element simulations of the plate-impact tests based on the model proposed by themselves and compared the numerical results to the experiments on tantalum. An increase of the spall strength was observed with higher shock intensities in the simulations.

Meanwhile, in the community of computational mechanics, a number of hydrocodes by which high-velocity impact calculations are performed have appeared in the last 30 years. These codes can be classified basically into three kinds: Lagrangian formulation (such as TOODY, EPIC and DYNA), Eulerian formulation (such as CSQ, CTH and MESA) and arbitrary Lagrangian–Eulerian (ALE) formulation (such as CAVEAT and SALE) [13]. Anderson [14] summarized the advantages and disadvantages of both Lagrangian and Eulerian descriptions. Johnson and Anderson [15] gave a review of the development of hydrocodes from a historical point of view. Benson [13] published a comprehensive survey which describes in detail the algorithms that are currently used in the modern hydrocodes. Lagrangian schemes work well at multimaterial interfaces, but suffer from their own troubles in problems with large deformations. The extreme case is that the simulation fails when a grid cell folds over itself. The way to avoid computational failure is to remesh the scheme where large distortion appears, which decreases the accuracy of the solution and the computational efficiency. Although the ALE method allows the mesh to move, the adaptive remeshing technologies are also required in regions with large deformations. Eulerian codes, however, allowing the boundaries to flow through a fixed mesh while computing the flow field on the fixed mesh, can solve arbitrarily large distortions. Coupled with an appropriate interface capturing algorithm, the Eulerian formulation will be quite attractive to high-velocity impact problems with large deformations.

Traditional Eulerian codes, mainly using first-order difference schemes and the von Neumann artificial viscosity, have low accuracy and cannot capture strong shock waves very well. In recent years, some researchers have applied the high-resolution techniques in the field of modern computational fluid dynamics (CFD) to the elastic–plastic flow problems. Benson [16] proposed an Eulerian formulation which is of the Lagrangian-plus-remap type and applied it to impact problems involving penetration and fracture. Cooper et al. [17] extended the Lagrangian-plus-remap method to simulate the dynamic void collapse in copper to study hot spot formation and jetting. Udaykumar et al. [18] presented a technique for the numerical simulation of high-speed multimaterial impact. The computations were performed on a fixed Cartesian mesh by casting the equations governing material deformation in Eulerian conservation law form. A high-order accurate ENO scheme was adopted along with the interface tracking technique to evolve sharp boundaries. After that, Tran and Udaykumar [19] updated their methodology by substituting the hybrid particle level set method [20] for the interface tracking algorithm adopted in the previous paper. These efforts have not only improved the accuracy of the numerical results for the problems of elastic–plastic flow, but also shown a remarkable direction in the field of computational mechanics.

The level set method (LSM), as mentioned above, was first devised by Osher and Sethian [21] as a simple and versatile method for tracking the boundary of a moving interface. The principle of this method is to describe an interface by the zero level set of a smooth function, called the level set function, and to update this function with Hamilton–Jacobi equations knowing the speed of the interface. The development and the application of the LSM were discussed in a review article by Osher and Fedkiw [22]. Enright et al. [20] improved the mass conservation properties of the LSM by using Lagrangian marker particles to rebuild the level set in regions which were being under-resolved. The method proposed by them is called the hybrid particle LSM.

A few years ago, Stolarska et al. [23] introduced the LSM to the field of the numerical simulation of crack propagation in solids. They coupled the LSM with the extended finite element method (FEM) and obtained satisfying results. Lately, Duflot [24] has reviewed in detail the update techniques of the crack representation by level set functions and proposed several new techniques. In these methods, a crack is an open curve (an open surface in three dimensions) that grows from its tip (its front in three dimensions), and two level set functions are necessary to represent a crack. When it comes to the spall fracture, however, this kind of crack representation becomes inappropriate. The crack appearing in the spall fracture opens wide with the motion of the scab. If the crack is modeled as a curve, the scab will never separate. In addition, it is unreasonable to put an initial crack in the target plate before simulation. In the present paper, a novel representation of crack by level set is proposed. The crack is represented as a two-dimensional narrow region which can be in any shape according to the specific problem. Free surfaces are created due to the crack.

The space–time conservation element and solution element (CE/SE) method, originally devised by Chang [25], [26], is a novel numerical scheme for solving hyperbolic conservation laws. It has several attractive features: (a) being mathematically simple; (b) a unified treatment of both space and time and enforcement of flux conservation in both space and time; (c) very little or almost no numerical dissipation; (d) the lack of directional splitting for flows in multiple spatial dimensions, resulting in a truly multidimensional scheme. The features mentioned above make the method substantially different from traditional well-established methods such as the finite difference and the finite volume methods. The CE/SE method has been generalized to high-order accuracy [27] and has reached a great success in the field of CFD [28], [29], [30], [31], [32], [33], [34]. Loh et al. [28] tested the CE/SE scheme for several problems ranging from linear acoustic waves to strongly nonlinear phenomena, with special emphasis on mixing-layer instability, and obtained satisfying numerical results. Wang et al. [32] simulated the detonation propagations with detailed chemical reaction models in a second-order accuracy CE/SE method. Recently, Wang et al. [35] have proposed an improved CE/SE scheme and extended it to the community of high-speed impact dynamics in solids. They have simulated two classical impact problems, Taylor bar impact and penetration. The comparisons with the experimental data and the results from other literature have proved the reliability of the computational technique developed by them.

The present paper extends the methodology proposed by Wang et al. [35] to simulate the spall fracture by introducing the idea of “element erosion” from the FEM. We apply the improved CE/SE scheme [35] to the high-velocity impact problems containing elastic–plastic flows, high strain rates and spall fractures in an axisymmetric framework. The Steinberg–Guinan (SG) constitutive model [36] and the Mie-Gruneisen (MG) equation of state (EOS) [37] are used to describe the target plate. The hybrid particle LSM [20] is adopted for capturing material interfaces and the physical boundary conditions are applied at the interfaces. We investigate the free boundary conditions and the contact boundary conditions in detail on the stress, velocity, density and internal energy. The VG model [5] is introduced to the codes to monitor the void volume fraction (VVF) of each element. The element is deleted, when its VVF reaches a critical value, to express the formation and propagation of a crack. Simulations of plate-impact experiments performed by Russian scientists [38] are presented. Our numerical results are in good agreement with the experimental observations described in the Russian spall database, which proves the feasibility and reliability of our computational technique.

Within the authors’ scope of knowledge, it is the first time to apply the LSM to describe the spall fracture. Since the Eulerian formulation can solve arbitrarily large distortions and the LSM has been successfully used for capturing interfaces undergoing extreme topological changes, the methodology presented here has advantages in simulating dynamic fracture such as spall in ductile materials caused by high-velocity or hypervelocity impact.