Abstract
We present a meshfree approach to model dynamic fracture in thin structures. Material failure is modeled based on a stress-based criterion and viscoplastic is used to describe the material behavior in the bulk material. Material fracture is simply modeled by breaking bonds between neighboring particles. The method is applied to fracture of cylindrical thin structures under explosive loading. The loading is modelled by a pressure-time history. Comparisons between the computational results and experimental data illustrate the validity and robustness of the proposed method.
1. Introduction
Modeling dynamic fracture of thin-walled structures remains a challenging task in computational mechanics. Such applications are of major importance in civil engineering, aeronautical engineering, aerospace engineering, and mechanical engineering. Thin structures are often modelled by shell theory. When shear effects can be neglected, the Kirchhoff-Love (KL) condition requires continuity of the underlying discretization. Effective formulations exploiting the higher order continuity of the meshfree methods have been exploited by [1, 2] for the first time; see also the contribution by [3–5] in the context of Isogeometric Analysis. The contribution in [1, 2] includes also fracture of the thin structure under dynamic loading. In [6], a fully coupled fluid-structure interaction method for fracturing thin structure has been proposed. On the other hand, there are many applications where shear effects need to be accounted for. There are numerous finite element formulations based on Mindlin-Reissner theory; see, for example, the manuscripts by [7–21]. On the other hand, thin-walled structures can also be modelled by three-dimensional continuum modeling.
Meshfree methods [22–34] are promising alternatives to finite element methods for applications involving large deformations, fracture and fluid-structure interaction. Meshfree method commonly exploits Lagrangian formulation [35–43] for both fluid and structure but does not suffer drawbacks of Lagrangian based finite element methods. They can handle large deformation and also fracture in a natural manner [44–51]. An excellent review of meshfree methods is given in [52, 53].
While there are many meshfree formulations for continua, there exist far less meshfree formulations for structures. A meshfree thin shell formulation based on KL theory has been presented by [54]. They employed the element-free Galerkin (EFG) method. Krysl et al. did not consider fracture. The first meshfree thin shell method for fracture was proposed by Rabczuk et al. [2]. This method was developed for linear and nonlinear materials and applied to static and dynamic problems. Moreover, it was developed for a finite-strain framework. The crack was modelled by partition-of-unity enrichment. Subsequently, [55–57] have simplified the treatment of cracks in thin shell and exploited an extrinsic basis [1]. However, the ill-posedness due to the enrichment scheme has not been addressed. Donning and Liu [58] used the reproducing kernel particle method (RKPM) [23] to model beam and plate elements. This approach has later been extended by [59]. Membranes and shells based on higher order approximations have been implemented by [60]. They report the higher order functions that are necessary to avoid locking. Wang and Chen [61] presented a locking-free meshfree formulation and showed that the Kirchhoff mode in the Mindlin plate can be reproduced by second-order approximations. Yagawa and Miyamura [62] combined a mixed approach within a free mesh framework. They used discrete Kirchhoff theory. Other approaches for thin structures based on 3D continuum approaches have been developed, for example, by the group of Li et al. [63].
As we want to keep the methodology applicable for thin structures with and without shear deformations, we will use a 3D continuum formulation that can be used also for real 3D simulations. The method will be applied to fracturing cylindrical structures under explosive loading. As the interaction between the solid and the fluid is highly complex [64–66], we will use a pressure-time history for the loading. The loading will be estimated by explosion simulations in rigid cylinders. Material failure of the solid is realized by simply breaking bonds between neighboring particles once a certain threshold is exceeded.
The paper is structured as follows. We first present the meshfree method, then briefly the constitutive model, and finally the fracture criterion. Then, we apply the method to fracturing thin cylinders. The numerical simulations will be compared to experimental data. The paper will be summarized in the last section.
2. Weak Form and Meshfree Discretization
Let us consider the linear momentum equation in the weak form. It is given in a total Lagrangian description of motion by finding the displacement field for all such that the variation in the energy is zero as with where and are the internal and external energy, respectively, denotes the kinetic energy, are the material coordinates, is the domain in the initial configuration and is the boundary in the initial configuration, denotes the body force, is initial density, is displacement, is the first Piola-Kirchhoff stress tensor, is the applied traction, denotes spatial derivatives with respect to material coordinate, and superimposed dots denote material time derivatives.
The boundary conditions are given by with boundaries and , the index referring to traction boundaries and the index to displacement boundaries; is the normal vector to the traction boundary.
The discretization of the displacement field in the meshfree approximation is given by where denotes the meshfree shape functions and denotes the number of particles. The EFG method [22] is based on a moving least squares approximation. It can be shown that the shape functions are expressed by with the moment matrix where the matrix contains the polynomial basis . We employ linear basis polynomials to ensure that the method fulfills the patch test. The matrix contains usually rational weighting functions . For computational efficiency, usually weighting functions with compact support are chosen; is the dilation parameter determining the size of the domain of influence of a central particle. We use the quartic B-spline function: with .
The test and trial functions have the structure of (4). Substituting them into the weak formulation leads to the following equation: After some algebraic operations, the well-known equation of motion is obtained as with
3. Material Model
A viscoplastic constitutive model from Zhou et al. [67] is used to describe the material behavior in the bulk. The material model in the bulk is of utmost importance for reliable predictions of material and structural behavior. Thermal effects are incorporated into the constitutive model as they are crucial in dynamic fracture simulations of thin structures as presented in this paper. The rate equations can be expressed in terms of the Jaumann rate of the Kirchhoff stress by where is the fourth order elasticity tensor, is the symmetric part of the velocity gradient , is the Jaumann rate of the Kirchhoff stress, with being the antimetric part of the velocity gradient, is the thermal expansion coefficient, and is the second-order identity matrix. The viscoplastic overstress model here is based on von Mises as with with being the back stress which is assumed to be zero here. A power low governs the thermal viscoplastic flow by with In (15) and (16), is a reference strain rate, is the rate sensitivity parameter, is the yield stress, is the corresponding reference strain and is Young’s modulus, is the strain hardening exponent, is a reference temperature, and and are thermal softening parameters. The function is the stress-strain relation measured at quasistatic strain rate of at temperature . The equivalent plastic strain is defined as Softening in material due to temperature is accounted for by varying material parameters where and are Young’s modulus and the Poisson ratio at temperature .
We employ the tangent modulus approach to update the thermal viscoplastic constitutive model. The details of the scheme can be found in the excellent manuscript by [68].
A stress-based criterion is employed to model fracture. Connectivities between neighboring nodes are broken when the maximum principal tensile strength is three times the strength . More sophisticated models based, for example, on loss of material stability might be used in future research.
4. Results
We validate our computational approach by comparison with experimental results. Therefore, we compare our numerical results to the fracture experiments by [69] who studied the failure mechanism of thin aluminum cylinders under detonation loading. They found that the failure mechanism changes with the length of the longitudinal precracks. Classical test set-up of their experiments is illustrated in Figure 1. The length of the initial cracks varies from 1′′ to 3′′:(i)1′′ (short notch),(ii)2′′ (medium-length notch),(iii)3′′ (long notch). In the experiments, the cylinder is filled with a combustible gas consisting of oxygen and ethylene. The pressure after ignition varies between 80 kPa and 180 kPa. After the gas is ignited at the left end, a detonation wave travels through the cylinder. In the thin aluminum specimen, pressures close to the Chapman-Jouguet (CJ) limit were measured. Its velocity is between 2300 m/s and 2400 m/s and the pressure values in the fully recreated CJ state range from 2.6 MPa to 6.1 MPa (depending on the initial pressure).
As it is computationally expensive to carry out a fully coupled FSI simulation as, for example, done in [6], we first study the pure detonation in a rigid tube needed to estimate the pressure-time history conditions in the aluminum cylinder. Those simulations are also based on meshfree simulations where the fluid is modelled with the element-free Galerkin method. A single exothermic chemical reaction A A B with a progress variable corresponding to the mass fraction ratio between the partial density of the reactant and the total density , that is, , governs the detonation by a modified Euler equation that contains additional inhomogeneous conservation laws: The chemical reaction follows the Arrhenius law [70] with being obtained from the following equation of state: where the parameter denotes the heat release due to the chemical reaction per unit mass [70]. The volume burn model according to [71] is exploited in order to guarantee the correct propagation speed and state in chemical equilibrium at all points in the discretization. Figure 2 illustrates the pressure-time history exemplary for one point close to the precrack. The pressure at the entire interior surface is monitored during the course of the simulation and subsequently applied as pressure-time history to study the fracture behavior of thin aluminum cylinders. We are aware that effects due to the changing structure after deformation are not included in our model. These effects might have a significant influence in particularly on fatal failure as the results below indicate. As stated previously, a fully coupled FSI model might be inevitable in this case.
For the fracture simulations, we employ discretizations of different refinements ranging from 140,000 particles to 280,000 particles. The following material parameters are used for all simulations: GPa, kg/m3, , GPa, , , and .
4.1. Short Notch
Figure 3 shows the deformed cylinder for a short notch. The crack grows straight in longitudinal direction. At a certain point, it curves in a 45-degree angle and continues to propagate in circumferential direction before it is arrested. At the end of the simulation, the fracture encompasses approximately 75% of the circumference. Similar results were obtained in the experiments from Chao [69] suggesting that effects due to fluid-structure interaction are not relevant.
| (a) Simulation; |
(b) Experiment
5. Medium-Length Notch
Different failure pattern is observed for the medium-length notched cylinder. In the experiments, the fracture farther from the detonation wave propagated first. After initial straight crack propagation, the fracture branched in 2 circumferential cracks separating the specimen in 2 pieces; see Figure 4.
| (a) |
| (b) |
| (c) |
(d) Experiment
Our computational results also predict severe failure though we are not capable of capturing the complete separation of the cylinder into two pieces. The fractures grow in two directions simultaneously and quickly turn into radial directions. The crack opening is more pronounced as for the short notch though. However, the fracture does not encompass the circumference. We believe that FSI simulations might be able to capture the failure mechanism better.
6. Long Notch
When the notch length is increased to 3′′, the specimen was broken into three pieces (in the experiments). This behavior is indeed captured well by our numerical simulations. Fractures grow first in longitudinal direction again. Several branches evolve leading finally to the failure illustrated in Figures 5 and 6.
| (a) |
| (b) |
| (c) |
| (d) |
| (a) |
| (b) |
(c) Experiment
7. Conclusions
In conclusion, we presented a simple and robust meshfree method for dynamic fracture in thin structures. The method is based on a 3D continuum approach and a simple fracture criterion is used to nucleate or propagate a crack. A viscoplastic model is employed in the bulk.
The method was applied to dynamic fracture in thin aluminum cylinders. The experiments from [69] were exploited to validate our method. In those experiments, a shock wave initiated by an explosive mixture was propagating through the aluminum cylinders. Different notch lengths were studied in their experiments, ranging from 1′′ to 3′′. As it is computationally and extremely challenging to model the entire process, we first performed simulations of the explosion in rigid tubes. The pressures in the rigid tubes were monitored during the course of the simulation and in a second step applied as boundary conditions to the thin aluminum cylinder, the specimen of interest. In other words: our simulations neglect effects due to fluid-structure interaction. Nevertheless, for the 1′′ and 3′′ notched specimen, our numerical simulations were able to capture the appropriate failure mechanisms while for the 2′′, we observed discrepancies to the experimental results.
Future studies will therefore focus on the extension of our model to fluid-structure interaction. Moreover, we intend to incorporate more advanced fracture models as proposed, for example, by [39, 72–75].
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this article.