Investigation of the Fracture and Fragmentation of Implosively Driven Thin-Walled Cylindrical Shell: From Thermodynamic Analysis to CDEM Simulation

The scattering of fragments is a notable characteristic of the explosive detonation of a shelled charge. This study examines the fracture and fragmentation of the shell and the process by which natural fragments form under the strains of implosion. The analysis takes into account both the explosive’s energy output and the casing’s dynamic response. For this purpose, utilizing a thermochemical code as an alternative to the conventionally employed cylinder test, the Jones–Wilkins–Lee equation of state (JWL EOS) was calibrated within a range of relative specific volume up to 13. The detonation of the shelled charge was subsequently analyzed using the continuum–discontinuum element method (CDEM). Following this, the formation mechanisms and scattering characteristics of natural fragments were scrutinized. The analysis found that the shell predominantly experiences shear failure with uniform evolution, displaying a “hysteresis effect” and two mutation stages in the evolution of tensile failure. Within the JWL EOS’s calibrated range, the representation of fragment displacement and velocity improved by 47.97% and 5.30%, respectively. This study provides valuable guidance for designing the power field of warheads and assessing their destructive power.


Introduction
Fragments generated by shelled charge explosions exhibit certain attributes, namely, rapid speed, minute size, large quantity and high destructive force.To forecast uncontrolled warhead detonations and aid in the design of controlled warheads, it is imperative to research the fracture laws and scattering characteristics of shells [1][2][3].
The velocity and mass-quantity distribution of natural fragments, which are principal measures of destructive potency, are intrinsically linked to the energy output of explosives and the dynamic mechanical characteristics of the shell [4][5][6].The equation of state for detonation products delineates the former, while a shell's material and structural parameters set the latter.The interaction of these factors determines how natural fragments form and scatter.
The equation of state, which depicts the relationships between the parameters of detonation products under the thermodynamic equilibrium state, is a crucial representation of an explosive's work capacity.Numerous forms of the equation of state have been developed.The semi-empirical equations of state, such as JWL EOS, which do not consider product composition change, are simple yet dependably accurate within certain temperature and pressure ranges, making them widely used in finite element calculation programs [7].However, their application is constrained because of over-simplification of the detonation problem from a phenomenological perspective.Thermochemical equations

Calculation of Chemical Reaction Perspective
The work capacity of explosives is substantially influenced by the thermodynamic parameters of detonation products such as pressure and temperature.Therefore, a comprehensive evaluation of explosive performance requires further procurement of CJ detonation parameters such as detonation velocity, pressure, temperature, and entropy.Using the DLCHEQ, the CJ point state parameters of the explosives were computed.Thereafter, the isentropic pressure-volume (P-V) relationship of the detonation products was calculated by leveraging the values of the entropy and detonation pressure.

Thermochemical Code
DLCHEQ is a thermochemical code developed based on the fundamental principles of thermodynamics and thermochemistry [37,38].It demonstrates reliable precision in describing the equation of state of detonation products over a wide spectrum of temperature and pressure.This capability allows for the comprehensive characterization of detonation products spanning from tens of GPa to atmospheric pressure and from thousands of K to room temperature.DLCHEQ can be divided into three core sections: the unreacted explosive equation of state, the EOS of product, and the model for reaction rate.According to different pressure ranges, the MCR equation (applicable from tens of GPa to hundreds of MPa), the second-order Virial equation (applicable from hundreds of MPa to atmospheric pressure), and the ideal gas equation (applicable from below atmospheric pressure) are applied in DLCHEQ.The EOS of detonation products encompasses sixteen distinct components, including CO 2 , O 2 , H 2 , H 2 O, CH 4 , NH 3 , NO, N 2 , CO, N, O, N 2 O, NO 2 , C, Al, and Al 2 O 3 .Through the input parameters of the explosive's molecular formula, the molecular formula of each detonation product component, and the explosive's initial density, DLCHEQ effectively generates the essential thermodynamic quantities of the detonation products, including volume, internal energy, composition, and enthalpy.Depending on the pressure, different equations of state are chosen to describe different gas components.After the temperature, pressure, and constituent components are identified, all the thermodynamic parameters of the detonation product system can be calculated through mass weighting.The sequence of these calculations is encapsulated in Figure 1.In assessing the states of chemical equilibrium and related calculations, the Gibbs free energy reaches its minimum value amid the chemical equilibrium condition.The chemical equilibrium equations are introduced as follows: In assessing the states of chemical equilibrium and related calculations, the Gibbs free energy reaches its minimum value amid the chemical equilibrium condition.The chemical equilibrium equations are introduced as follows: where α ik denotes the number of atoms of the kth element in the molecular formula of the ith component, λ k denotes the chemical potential of the ith element, n i and n0 denote the numbers of moles before and after the iteration of the ith kind, respectively, and n and n 0 denote the total number of moles before and after the iteration of the gas phase product, respectively.The calculation formulas are as follows:n denotes the free energy of the ith kind of gas phase product and G i (n 0 ) denotes the free energy of the ith kind of solid phase product.Through DLCHEQ, the isentropic expansion process of detonation products from the CJ state to the state of relative specific volume v = 13 were computed, along with the compositional evolution of detonation products and the isentropic expansion line.
This article focuses on a mixed explosive comprising 70% HMX (C 4 H 8 N 8 O 8 ) and 30% RDX (C 3 H 6 N 6 O 6 ) as the subject of calculation.The enthalpies of formation for the individual explosives are 0.255 kJ/g and 0.277 kJ/g, respectively, while their respective molecular weights are 296.155and 222.116.

Compositional Evolution
Spanning the range of relative specific volume from 1 to 13, the compositional evolution of detonation products is illustrated in Figure 2. Due to the order of magnitude difference of content, constituent components were represented by solid and dotted lines.As can be seen, the system is predominantly composed of N 2 , H 2 O, CO 2 , and CO, with lesser quantities of NH 3 , CH 4 , and H 2 .The concentration of all component parts fluctuates noticeably, with the exception of N 2 .The H 2 O content diminishes consistently during expansion, while C and CO 2 react to yield CO. When the C (graphite) concentration diminishes to zero, the increase in CO transitions to a decrease simultaneously with an increase in CO 2 .This suggests that the reaction experiences an initial exothermic process followed by heat absorption.noticeably, with the exception of N2.The H2O content diminishes consistently during expansion, while C and CO2 react to yield CO. When the C (graphite) concentration diminishes to zero, the increase in CO transitions to a decrease simultaneously with an increase in CO2.This suggests that the reaction experiences an initial exothermic process followed by heat absorption.
The constituent components undergo an overall evolution characterized by an initial rapid progression followed by deceleration.The evolution process of the constituent components is divided into two distinct parts demarcated by the measurable limits of the 25 mm cylinder test ( v = 7).The first part reveals the occurrence of the primary reaction, while the reaction in the second part continues at a reduced rate.Comparing the end points of part 2 with part 1, the percentage of CO is 62.48% of the corresponding value, while CO2 amounts to 111.20% and CH4 amounts to 122.97% of their respective values in part 1.These findings illustrate the persistence of the reaction between detonation products beyond the measurable range of the cylinder test, exerting an influence on the compositional evolution.The work capacity of explosives is influenced by the compositional evolution of the detonation products.However, in the process of formulating the empirical equation of state, the constituent compositions are usually regarded as unchanged.The range of meas- The constituent components undergo an overall evolution characterized by an initial rapid progression followed by deceleration.The evolution process of the constituent components is divided into two distinct parts demarcated by the measurable limits of the 25 mm cylinder test (v = 7).The first part reveals the occurrence of the primary reaction, while the reaction in the second part continues at a reduced rate.Comparing the end points of part 2 with part 1, the percentage of CO is 62.48% of the corresponding value, while CO 2 amounts to 111.20% and CH4 amounts to 122.97% of their respective values in part 1.These findings illustrate the persistence of the reaction between detonation products beyond the measurable range of the cylinder test, exerting an influence on the compositional evolution.
The work capacity of explosives is influenced by the compositional evolution of the detonation products.However, in the process of formulating the empirical equation of state, the constituent compositions are usually regarded as unchanged.The range of measurement in the cylinder test is restricted by the tube breaks.As shown in Figure 2, the constituent components continue to change substantially during the evolution process after v = 7, directly affecting the work capacity of the explosives due to endothermic and exothermic reactions.Hence, opting for the thermochemical code is a better choice for calibrating the JWL EOS as opposed to relying on the cylinder test.This approach affords extended ranges of temperature and pressure for the expansion data of detonation products, leading to a more precise depiction of the work capacity of explosives.

Isentropic Expansion Curve
The working range of detonation products is assumed as isentropic expansion.According to the computed entropy and detonation pressure, the isentropic P-V relationship of detonation products can be computed by DLCHEQ.
The shock compression state of the detonation products satisfies the equation of state; then, the Hugoniot relation is introduced as follows: where ρ denotes the particle number density.The pressure P and internal energy E of the system are obtained from the total free energy F, while E, P, and V are indicated by the subscripts 0 and H for the molar internal energy, pressure, and molar volume of the initial and shock-compressed final states of the system, respectively.Because the detonation velocity of the CJ state is at the minimum, the CJ detonation velocity can be solved by the parabolic minimum method for the three points on the Hugoniot curve.The calculation formula is introduced as follows: The main constituent components of the test charge in this paper are 70% HMX and 30% RDX.The CJ detonation parameters calculated by the above equation are shown as follows: detonation velocity D CJ = 9.1772 km/s, specific volume v CJ = 0.4164 cm 3 /g, detonation temperature T CJ = 3818.2K, detonation pressure P CJ = 35.286MPa, internal energy E CJ = 2.266 kJ/g, entropy S = 7.3016 J/g/K.In this way, the isentropic expansion curve of detonation products can be obtained.

Calibration of JWL Equation of State
The coefficients were calibrated using JWL EOS to fit the isentropic curve computed through DLCHEQ.The two-dimensional fluid dynamics software DYNA-2D was employed to simulate the cylinder test to obtain the expansion curve of the external wall.The reliability of the coefficients derived from the calibration was validated through comparison with the results of the explosive cylinder test in the range of v < 7.

Curve Fitting
The fitting function w+1) is the form of the isentropic JWL EOS, where P denotes the instantaneous pressure of the detonation product.Only three of the six coefficients A, B, C, R 1 , R 2 , and ω to be calibrated are independent.A genetic algorithm was used to fit R 1 , R 2 , and ω, while A, B, and C were solved simultaneously according to three compatibility equations.The compatibility equation is established by the CJ condition, Hugoniot relation, and −(∂p s /∂V)v CJ = ρ 0 D 2 of the ideal detonation state, which are introduced as follows: where the initial internal energy E 0 =10 GPa, the range of values of R 1 , R 2 , and ω is R 1 = 4-5, R 2 = 1-2, and ω = 0.2-0.4,respectively.The fitting results are shown in Figure 3.The maximum error on the isentropic curve ∆P = 0.6 GPa.A closer inspection of the curve's intermediate section reveals a satisfactory fitting, suggesting that the genetic algorithm applied in this study accurately fits the product's isentropic expansion curve.Coefficients of the JWL EOS obtained from the calibration are shown in Table 1.A 25 mm cylinder test was conducted in the explosion tower, utilizing a xenon lamp for lighting purposes.Data were procured from an ultra-high-speed scanning camera operating at a velocity of 3 mm/μs and subsequently calibrated with the assistance of PDV laser interferometry.The designated observational vantage point was positioned 120 mm distant from the cylinder's terminus.The sample, instrument and tower layout are shown in Figure 4.The test cylinder, composed entirely of TU1 oxygen-free pure copper (density of 8.93 g/cm 3 ), exhibited a charge density of 1.86 g/cm 3 .Coefficients of the JWL EOS obtained from the calibration are shown in Table 1.

Test Verification
A 25 mm cylinder test was conducted in the explosion tower, utilizing a xenon lamp for lighting purposes.Data were procured from an ultra-high-speed scanning camera operating at a velocity of 3 mm/µs and subsequently calibrated with the assistance of PDV laser interferometry.The designated observational vantage point was positioned 120 mm distant from the cylinder's terminus.The sample, instrument and tower layout are shown in Figure 4.The test cylinder, composed entirely of TU1 oxygen-free pure copper (density of 8.93 g/cm 3 ), exhibited a charge density of 1.86 g/cm 3 .DYNA-2D two-dimensional fluid dynamics software was utilized to simulate the cylinder test.Both the Johnson-Cook constitutive model and the Gruneisen equation of state [39] were applied for the copper tube.The simulation results are presented in Figure 5, illustrating the expansion curve of the cylinder wall.The JWL EOS calibrated by the isentropic curve shows significant consistency with the test measurements, which can effectively describe the expansion process of the external wall in the cylinder test.Thus, the reliability of the JWL EOS, as calibrated in this study is substantiated, asserting its capability to accurately characterize the work capacity of explosives.DYNA-2D two-dimensional fluid dynamics software was utilized to simulate the cylinder test.Both the Johnson-Cook constitutive model and the Gruneisen equation of state [39] were applied for the copper tube.The simulation results are presented in Figure 5, illustrating the expansion curve of the cylinder wall.The JWL EOS calibrated by the isentropic curve shows significant consistency with the test measurements, which can effectively describe the expansion process of the external wall in the cylinder test.Thus, the reliability of the JWL EOS, as calibrated in this study is substantiated, asserting its capability to accurately characterize the work capacity of explosives.DYNA-2D two-dimensional fluid dynamics software was utilized to simulate the cylinder test.Both the Johnson-Cook constitutive model and the Gruneisen equation of state [39] were applied for the copper tube.The simulation results are presented in Figure 5, illustrating the expansion curve of the cylinder wall.The JWL EOS calibrated by the isentropic curve shows significant consistency with the test measurements, which can effectively describe the expansion process of the external wall in the cylinder test.Thus, the reliability of the JWL EOS, as calibrated in this study is substantiated, asserting its capability to accurately characterize the work capacity of explosives.The accuracy of the JWL EOS was verified over the range of relatively specific volumes measured in the cylinder test.The 25 mm cylinder test can describe the work capacity of the explosive in the range of v < 7, while the explosive continues to exert a propulsive acceleration effect on fragments in v > 7. The thermochemical code DLCHEQ computes the isentropic expansion curve over a wider temperature-pressure range.This process takes into account both the composition of the explosive and the chemical reaction process.Consequently, the JWL EOS can characterize the explosive's work capacity over a longer time interval.

Thin-Walled Cylindrical Shells Driven by Detonation Products
Driven by detonation products, the thin-walled shell fractures and natural fragments are subsequently produced.Experiments and corresponding numerical simulation were conducted to investigate the fracture characteristics of shells and acceleration of fragments.The velocities of the fragments were measured in the experiment, and the JWL EOS of detonation products obtained from the previous section was applied to the numerical simulation utilizing CDEM.Simulation results of the initial velocities and mass-quantity distribution of the fragments were validated with experimental results and Mott theory.

Experiment of Shelled Charge Explosion
A shelled charge explosion test was conducted to investigate the velocity and mass distribution of fragments.A Q235 steel shell with the following specifications was used in the experiment: length L = 285 mm, inner diameter D I = 84.2mm, outer diameter D O = 88.8 mm, shell wall thickness δ = 2.3 mm, shell mass M S = 1.36 kg, charge mass M E = 3 kg.Fragment velocity was measured electrically with a comb-shaped target.There is a short-circuit in the velocity signal output when the fragment does not reach the target; conversely, the signal output is an open circuit when the fragment arrives at the target.This system offers accurate measurement of fragment velocities at various spatial locations, thereby mitigating the constraints of the traditional single-point through-break type velocity target.The impact of explosion and the scattering effect of fragments are illustrated in Figure 6.The energy release framework in a shelled charge explosion encompasses shockwaves, fragment dissemination, high-temperature fireballs, etc.
Materials 2023, 16, x FOR PEER REVIEW 9 of 24 ity of the explosive in the range of v < 7, while the explosive continues to exert a propul- sive acceleration effect on fragments in v > 7. The thermochemical code DLCHEQ com- putes the isentropic expansion curve over a wider temperature-pressure range.This process takes into account both the composition of the explosive and the chemical reaction process.Consequently, the JWL EOS can characterize the explosive's work capacity over a longer time interval.

Thin-Walled Cylindrical Shells Driven by Detonation Products
Driven by detonation products, the thin-walled shell fractures and natural fragments are subsequently produced.Experiments and corresponding numerical simulation were conducted to investigate the fracture characteristics of shells and acceleration of fragments.The velocities of the fragments were measured in the experiment, and the JWL EOS of detonation products obtained from the previous section was applied to the numerical simulation utilizing CDEM.Simulation results of the initial velocities and mass-quantity distribution of the fragments were validated with experimental results and Mott theory.

Experiment of Shelled Charge Explosion
A shelled charge explosion test was conducted to investigate the velocity and mass distribution of fragments.A Q235 steel shell with the following specifications was used in the experiment: length L = 285 mm, inner diameter DI = 84.2mm, outer diameter DO = 88.8 mm, shell wall thickness δ = 2.3 mm, shell mass MS = 1.36 kg, charge mass ME = 3 kg.Fragment velocity was measured electrically with a comb-shaped target.There is a shortcircuit in the velocity signal output when the fragment does not reach the target; conversely, the signal output is an open circuit when the fragment arrives at the target.This system offers accurate measurement of fragment velocities at various spatial locations, thereby mitigating the constraints of the traditional single-point through-break type velocity target.The impact of explosion and the scattering effect of fragments are illustrated in Figure 6.The energy release framework in a shelled charge explosion encompasses shockwaves, fragment dissemination, high-temperature fireballs, etc.The velocity tests were executed at distances of 4 m and 5.2 m, forming three distinct test groups, with the collected data illustrated in Table 2. -0   is utilized for cal- culating the initial velocity of the fragment, where Vx denotes the velocity of the fragment when it is scattered by x m, V0 denotes the initial speed of the fragment, and a denotes the coefficient to be fitted.Substituting the data in Table 2 gives V0 = 2.22 km/s.The comb target, being a type of on-off velocity measuring system, measures the velocity of the fragment that first impacts the target surface.Thus, under the static explosion condition of this experiment, the maximum velocity of the fragment is 2.22 km/s.The ratio of charge to shell mass for the test, ME/MS = 2.21, is typical of thin-walled shells, implying a high The velocity tests were executed at distances of 4 m and 5.2 m, forming three distinct test groups, with the collected data illustrated in Table 2. V x = V o × e −ax is utilized for calculating the initial velocity of the fragment, where V x denotes the velocity of the fragment when it is scattered by x m, V 0 denotes the initial speed of the fragment, and a denotes the coefficient to be fitted.Substituting the data in Table 2 gives V 0 = 2.22 km/s.The comb target, being a type of on-off velocity measuring system, measures the velocity of the fragment that first impacts the target surface.Thus, under the static explosion condition of this experiment, the maximum velocity of the fragment is 2.22 km/s.The ratio of charge to shell mass for the test, M E /M S = 2.21, is typical of thin-walled shells, implying a high number of lighter, smaller fragments produced by the explosion.Considering the limitation of the recovery method, only a few parts of fragments were able to be recovered, with each weighing less than 3 g.These fragments account for a total mass less than that of the original shell, rendering it insubstantial as a comparative value.In the midst of detonation product expansion, the metallic shell undergoes three stages, respectively, deformation, micro-crack initiation, and expansion, resulting in formation of natural fragments.To accurately portray the process of crack initiation and propagation while maintaining the consistency of momentum, energy, and mass, the continuum-discontinuum element method was employed in this study to investigate the fracture and fragmentation characteristics of the metal shell subjected to implosion loading.

Basic Concept
The continuum-discontinuum element method (CDEM) can be defined as a dynamic explicit solution algorithm that functions within a Lagrangian framework.Strict control equations are established via the Lagrangian energy system, which, when combined with the dynamic relaxation method, allows for explicit iterative solutions.This enables the attainment of a unified depiction of continuum-discontinuum scenarios.In order to analyze material failure progression, the model takes into account fractures that occur both inside and on the boundaries of the block.The algorithm simulates the entire course, from continuous deformation of the material to eventual fracture and subsequent movement.This method integrates the benefits of both continuous and discontinuous calculation: the former implements techniques such as finite element, finite volume, and spring element, while the latter utilizes the discontinuous element method.
The numerical model in CDEM consists of two parts: block and interface.The block is composed of one or more finite element cells, which are used to characterize the continuous characteristics of the material such as elasticity, plasticity, and damage.The common boundary between the two blocks is the interface, which is used to characterize the discontinuous characteristics of the material such as fracture, slip, and collision.Within the CDEM framework, the interface includes two concepts: real interface and virtual interface.The real interface is used to characterize the real discontinuous surface, including the interface between constituent components.In contrast, the virtual interface serves dual functions.First, it connects two blocks to transmit mechanical information; second, it provides a potential channel for the propagation of explicit cracks (i.e., cracks can propagate along any virtual interface).The schematic diagram of the numerical model in CDEM is shown in Figure 7.The schematic model consists of seven blocks, of which one block is composed of two triangular elements and the remaining six blocks are composed of one triangular element each.In addition, the red line in Figure 7c is the real interface, while the black line is the virtual interface.

Explosion Model
The model applies the standard form of the JWL EOS to delineate the expansion process of detonation products.The calculation formula is introduced as follows:

Explosion Model
The model applies the standard form of the JWL EOS to delineate the expansion process of detonation products.The calculation formula is introduced as follows: where E 0 denotes the specific internal energy of detonation products at the initial moment (J/m 3 ).The values of A, B, R 1 , R 2 , and ω are shown in Table 1.
It is assumed that the ignition time of an explosive (including several units) is t 0 , the coordinates of the ignition point are (x 0 , y 0 , z 0 ), the distance from the centre of a unit in the explosive to the ignition point is d, and the detonation velocity of the explosive is D, making the ignition time of the unit t 1 = d/D + t 0 .When the explosion time t > t 1 , the explosion pressure is calculated according to Equation ( 6), where P r denotes the real explosion pressure; here, P(v, E) denotes the equation of state of detonation products based on the JWL model and ξ denotes the energy release rate, which can be obtained by Equation (7), where V e denotes the initial volume of the unit and A e-max denotes the maximum area of the unit.
Furthermore, this study presents a model for detonation product leakage, with the aim of depicting the decrease in pressure following shell break.Considering the inevitable dissipation of detonation gas post-explosion, an exponential decay function serves as an effective tool to articulate the decline in detonation gas pressure resulting from such dissipation.The calculation formula is introduced as follows: where p denotes the explosion pressure calculated according to the formula of the explosion source, t c denotes the characteristic time (when t = t c , the pressure becomes 36.8% of the original pressure.),n denotes the characteristic index, and t c is closely related to the geometry of the shell and the charge, which is chosen to be t c = 90 µs.

Constitutive Model
Based on the shell size and charge type in the explosion test, the numerical model is shown in Figure 8.It utilizes a tetrahedral grid of 2 mm size, with a total grid number of 3.6 × 10 5 .The test charge is represented by a blue cylinder with radius R I = 42.1 mm, while the grey cylindrical shell with radius R O = 44.4mm and wall thickness δ = 2.3 mm symbolizes the metal shell.These two components make uniform contact, with the length of the shell being L = 285 mm.
Five evenly spaced monitoring markers are positioned on the shell's surface to capture the information of fragments displacement and velocity at the designated location.The layout of these points is shown in Figure 8. Specifically, point P3 serves as the central point of the shell busbar, while points P1 and P2 are symmetric to points P5 and P4, respectively, relative to the radial center section of the shell.
In order to accurately simulate the fracture process of the shell, the bilinear fracture model is used to describe the fracture evolution process between metal elements.First, the elastic contact force is calculated by Equation (9), where F n and F s denote the normal and tangential contact forces, respectively, K n and K s denote the normal and tangential contact stiffness, respectively, and ∆d n and ∆d s denote the relative displacement increments in the normal and tangential directions, respectively.
Materials 2023, 16, x FOR PEER REVIEW 12 of 24 the grey cylindrical shell with radius RO = 44.4mm and wall thickness δ = 2.3 mm symbolizes the metal shell.These two components make uniform contact, with the length of the shell being L = 285 mm.Five evenly spaced monitoring markers are positioned on the shell's surface to capture the information of fragments displacement and velocity at the designated location.The layout of these points is shown in Figure 8. Specifically, point P3 serves as the central point of the shell busbar, while points P1 and P2 are symmetric to points P5 and P4, respectively, relative the radial center section of the shell.In order to accurately simulate the fracture process of the shell, the bilinear fracture model is used to describe the fracture evolution process between metal elements.First, the elastic contact force is calculated by Equation ( 9), where Fn and Fs denote the normal and tangential contact forces, respectively, Kn and Ks denote the normal and tangential contact stiffness, respectively, and Δdn and Δds denote the relative displacement increments in the normal and tangential directions, respectively.
) ( ) ( ) ( ) In order to calculate the progressive failure process of the material, the Mohr-Coulomb criterion and the maximum tensile stress criterion are introduced as follows: (1) tensile failure (2) shear failure where T denotes the tensile strength of the current time step, C denotes the cohesive force of the current time step, f denotes the internal friction angle, and Ac denotes the contact area.The structural softening model of tensile failure and shear failure is shown in Figure 9.In order to calculate the progressive failure process of the material, the Mohr-Coulomb criterion and the maximum tensile stress criterion are introduced as follows: where T denotes the tensile strength of the current time step, C denotes the cohesive force of the current time step, f denotes the internal friction angle, and A c denotes the contact area.
The structural softening model of tensile failure and shear failure is shown in Figure 9.In order to characterize the effect of the strain rate on the mechanical parameters of the metal, the relationship between the cohesion and tensile strength of the interface and the tangential strain rate and normal strain rate of the interface is introduced as follows in CDEM.
where C0 and T0 denote the static cohesion and tensile strength, respectively, γ and ε denote the tangential and normal strain rates on the interface, respectively, and α denotes the strain rate coefficient.
The mechanical parameters of the Q235 steel used in the numerical simulation are shown in Table 3.In order to characterize the effect of the strain rate on the mechanical parameters of the metal, the relationship between the cohesion and tensile strength of the interface and the tangential strain rate and normal strain rate of the interface is introduced as follows in CDEM.
where C 0 and T 0 denote the static cohesion and tensile strength, respectively, γ and ε denote the tangential and normal strain rates on the interface, respectively, and α denotes the strain rate coefficient.
The mechanical parameters of the Q235 steel used in the numerical simulation are shown in Table 3.

Verification of Numerical Results
The numerical simulation spanned a duration from 0 to 100 µs.The mass-quantity distribution of the natural fragments was counted.The resulting displacement and velocity characteristics of the fragments were recorded by marker points on the surface of the shell and the fracture process and damage/fracture type of the shell were recorded by cell monitoring.

Initial Velocity of Fragments
The velocity cloud diagram of fragments at 100 µs is shown in Figure 10.The resulting velocity distribution of the majority of these fragments falls within the range of 1.8-2.1 km/s.The maximum, average, and minimum value of the resulting velocity are V Max = 2.118 km/s, V average = 2.019 km/s, and V Min = 1.902 km/s, respectively.Based on the data of the velocimetric target in Table 2, the initial velocity of the fragment was calculated as V0 = 2.22 km/s.The errors of VMax, Vaverage, and VMin were 4.59%, 9.04%, and 14.32%, respectively.These numerical simulation results can effectively reflect the real situation of fragment scattering velocity in the test.

Mass-Quantity Distribution
Existing fragment recovery technology faces significant limitations, making comprehensive recovery of all fragments following shelled charge explosion a challenging task.The Mott distribution theory has been extensively applied to assess the mass-quantity distribution of natural fragments.Furthermore, preceding research findings have corroborated the accuracy of the Mott distribution.This paper aims to examine the accuracy of numerical simulations, which is achieved by undertaking a comparative analysis between the numerical simulation outcomes and the established Mott theory.
The shell examined in this study represents a quintessential thin-walled structure.It manifests fractures in two dimensions, subsequently leading to the formation of natural fragments.The mass-quantity distribution of the fragments is introduced as follows: where N(mf) denotes the number of fragments with a mass greater than mf and 2μ denotes the average mass of broken fragments.Incorporating the principles of Mott Theory, Gur- Based on the data of the velocimetric target in Table 2, the initial velocity of the fragment was calculated as V 0 = 2.22 km/s.The errors of V Max , V average , and V Min were 4.59%, 9.04%, and 14.32%, respectively.These numerical simulation results can effectively reflect the real situation of fragment scattering velocity in the test.

Mass-Quantity Distribution
Existing fragment recovery technology faces significant limitations, making comprehensive recovery of all fragments following shelled charge explosion a challenging task.The Mott distribution theory has been extensively applied to assess the mass-quantity distribution of natural fragments.Furthermore, preceding research findings have corroborated the accuracy of the Mott distribution.This paper aims to examine the accuracy of numerical simulations, which is achieved by undertaking a comparative analysis between the numerical simulation outcomes and the established Mott theory.
The shell examined in this study represents a quintessential thin-walled structure.It manifests fractures in two dimensions, subsequently leading to the formation of natural fragments.The mass-quantity distribution of the fragments is introduced as follows: where N(m f ) denotes the number of fragments with a mass greater than m f and 2µ denotes the average mass of broken fragments.Incorporating the principles of Mott Theory, Gurney and Sarmousakis proposed a µ formulation applicable for thin-walled cylindrical shells is introduced as follows: where K denotes the parameter related to the performance of the explosive in g 1/2 /cm 7/6 .Concerning the known empirical parameters of monomeric explosives and the ratio of the components used in this test, the value of K was taken as 0.293 and µ = 0.01.The resulting fragment mass-quantity distributions of the Mott theoretical calculation and the numerical simulation are shown in Figure 11.The fragments were divided into five weight categories for analysis, which were formed by considering the ratio of fragment count and total weight to that in various mass intervals.For each weight category, the counts of fragments computed by Mott's theoretical calculation and derived from numerical simulations were compared.The discrepancy between these values is presented in Table 4 As can be seen from Table 4, the error in four of these intervals is less than 20% when the simulation results are compared with those of the Mott theory.The interval between 0.01-0.0105g is an exception, with a maximum error of 30.86%, attributable to the constraints of the smallest grid dimension in the numerical model.The minimum detectable fragment mass is 0.0074 g; hence, predictions for a total number of fragments closer to this mass might exhibit inaccuracies.Mott theory and numerical simulation estimate the total number of fragments as 3.486 × 10 4 and 3.2827 × 10 4 , respectively, yielding a simulationerror margin of 5.83%.These numerical simulation outputs serve as appropriate reflections of the mass-quantity distribution of natural fragments.
Subsequent to this analysis, the implosion-induced fracture and fragmentation process of thin-walled shells can be accurately portrayed using continuum-discontinuum element method.The fragments were divided into five weight categories for analysis, which were formed by considering the ratio of fragment count and total weight to that in various mass intervals.For each weight category, the counts of fragments computed by Mott's theoretical calculation and derived from numerical simulations were compared.The discrepancy between these values is presented in Table 4: As can be seen from Table 4, the error in four of these intervals is less than 20% when the simulation results are compared with those of the Mott theory.The interval between 0.01-0.0105g is an exception, with a maximum error of 30.86%, attributable to the constraints of the smallest grid dimension in the numerical model.The minimum detectable fragment mass is 0.0074 g; hence, predictions for a total number of fragments closer to this mass might exhibit inaccuracies.Mott theory and numerical simulation estimate the total number of fragments as 3.486 × 10 4 and 3.2827 × 10 4 , respectively, yielding a simulation-error margin of 5.83%.These numerical simulation outputs serve as appropriate reflections of the mass-quantity distribution of natural fragments.
Subsequent to this analysis, the implosion-induced fracture and fragmentation process of thin-walled shells can be accurately portrayed using continuum-discontinuum element method.

Discussion
The detonation of a shelled charge progresses through three phases: initial shell deformation and fracture, acceleration of fragments driven by detonation products, and a deceleration phase attributable to air resistance.Both the initial deformation/fracture and the fragment acceleration stages transpire within a very brief period, on the order of microseconds, while the deceleration stage progresses over a longer time scale, on the order of milliseconds.The fracture and fragmentation process of the shell were analyzed based on the verified numerical simulation results in the CDEM program.

Fracture Mechanism of the Thin-Walled Cylindrical Shell
Under the load of implosion, the shell deforms, and breaks as pressure from the detonation products diminishes.The failure modes include both tensile and shear failure.Using the continuum-discontinuum element method, the virtual interface between cells slips and breaks under the influence of the detonation products.To investigate the mechanism of the shell's fracture evolution, the damage to the spring and the fracture area of the shell element under shear and tensile loads were quantified based on the results of numerical simulations.The spring area of shear damage/fracture is on the order of 10 −1 m 2 , while the spring area of tensile damage/fracture is on the order of 10 −4 m 2 .The difference between the two is 1000 times, indicating that the failure of thin-walled shells is dominated by shear failure.At 100 µs, the shear damage area is 0.59 m 2 , while the shear fracture area is 0.38 m 2 ; 64.41% of the blocks with shear damage evolve into shear fractures.The tensile damage area is 3.53 × 10 −4 m 2 , while the tensile fracture area is 3.28 × 10 −4 m 2 .The blocks with tensile damage evolve into tensile fractures with a ratio of 92.91%.
At 30 µs, the shear and tensile spring damage and fracture area account for 99.97%, 90.22%, 99.97%, and 97.44% of the time of 100 µs, respectively, indicating that the fracture of the shell mainly occurs in the time interval of 0-30 µs.Therefore, the shell fracture evolution study was carried out around this time interval.The shell fracture evolution process is plotted in Figure 12.

Fracture Mechanism of the Thin-Walled Cylindrical Shell
Under the load of implosion, the shell expands, deforms, and breaks as pressure from the detonation products diminishes.The failure modes include both tensile and shear failure.Using the continuum-discontinuum element method, the virtual interface between cells slips and breaks under the influence of the detonation products.To investigate the mechanism of the shell's fracture evolution, the damage to the spring and the fracture area of the shell element under shear and tensile loads were quantified based on the results of numerical simulations.The spring area of shear damage/fracture is on the order of 10 −1 m 2 , while the spring area of tensile damage/fracture is on the order of 10 −4 m 2 .The difference between the two is 1000 times, indicating that the failure of thin-walled shells is dominated by shear failure.At 100 μs, the shear damage area is 0.59 m 2 , while the shear fracture area is 0.38 m 2 ; 64.41% of the blocks with shear damage evolve into shear fractures.The tensile damage area is 3.53 × 10 −4 m 2 , while the tensile fracture area is 3.28 × 10 −4 m 2 .The blocks with tensile damage evolve into tensile fractures with a ratio of 92.91%.
At 30 μs, the shear and tensile spring damage and fracture area account for 99.97%, 90.22%, 99.97%, and 97.44% of the time of 100 μs, respectively, indicating that the fracture of the shell mainly occurs in the time interval of 0-30 μs.Therefore, the shell fracture evolution study was carried out around this time interval.The shell fracture evolution process is plotted in Figure 12.Shear damage initially manifests in the shell at 4 μs, while tensile damage becomes evident at 4.5 μs, resulting in a 0.5 μs temporal discrepancy, with shear damage being the predecessor.Both shear and tensile fractures initiate almost simultaneously, at 5.3 μs and 5.4 μs, respectively.According to the law of propagation of stress waves in solids and the tensile stress guidelines for shell fracture, cracks begin to appear when the stress system Shear damage initially manifests in the shell at 4 µs, while tensile damage becomes evident at 4.5 µs, resulting in a 0.5 µs temporal discrepancy, with shear damage being the predecessor.Both shear and tensile fractures initiate almost simultaneously, at 5.3 µs and 5.4 µs, respectively.According to the law of propagation of stress waves in solids and the tensile stress guidelines for shell fracture, cracks begin to appear when the stress system of the shell wall reaches a critical state.Radial cracks appear in the circumferential tensile stress area of the shell and propagate along the tensile stress area.Furthermore, they are unable to extend into the compressive stress area, resulting in crack formation on the shell's outer surface with subsequent propagation towards the inner surface.When the thickness of the compressive stress region existing in the inner wall of the shell is reduced to zero, the compressive stress of the inner wall of the shell is equal to the yield stress of the material.The crack propagates to the inner surface and the shell is completely broken.During this process, shear fractures occur in the direction of maximal shear stress when an exterior shell crack develops.The tensile fracture occurs perpendicularly to the principal stress and only after the crack extends to the inner surface.Consequently, due to the shell wall's thickness, the tensile failure occurs later than the shear failure.This temporal disparity constitutes the tensile failure's "hysteresis effect".
The shell mainly incurs block damage at 4-18 µs.The shear damage area is 0.588 m 2 and the tensile damage area is 3.52 × 10 −4 m 2 at 18 µs.At 18-30 µs, the shear damage area increases by 0.575% and the tensile damage area increases by 0.387%, indicating that the block damage area stops growing after the shell has completely broken.Most of the block fracture occurs simultaneously with block damage, though it continues to grow after 18 µs.In the three stages of 4-18 µs, 18-30 µs, and 30-100 µs, the area produced by shear fracture accounted for 71.046%, 21.895%, and 7.059% of the total value, respectively.Additionally, the area produced by tensile fracture accounted for 88.644%, 9.759%, and 1.597% of the total value, respectively.This indicates that after 18 µs the shear fracture area increased slowly until 100 µs, while the tensile fracture area ended before 30 µs.
The crack proliferation cloud map of the shell is shown in Figure 13.For convenience of observation, the ratio distance of the cloud map is set to 0. Red and green denote the shear and tensile failure, respectively.The left side shows the cloud map depicting the failure face piece, while the right side shows the cloud map illustrating the presence of surface cracks on the shell.As can be seen from the cloud map, the shell crack propagates along the axial direction to both ends.Detonation products reach the end cover of the shell at 16.3 µs, completing more than 99.5% of the damage area growth.At 32.4 µs, the shear fracture area of 92.94% and tensile fracture area of 98.41% are increased, with a shell expansion for the initial diameter of 1.8 times.
As seen in the selected part of the dashed line frame in Figure 12b, the growth of the tensile damage/fracture area undergoes two mutations, with the mutation occurrence intervals for damage and fracture overlapping.The first mutation occurs in the time range of 5.5-6.5 µs, corresponding to the failure cloud map in Figure 13a.At this stage, the detonation products swiftly reach the shell wall, initiating the formation of a crack along the circumference starting from the center of the shell.During this stage, the evolution time represents only 1% of the total time, while the growth area attributed to tensile damage accounts for 18.70% of the total area and the growth area associated with tensile fracture amounts to 24.91% of the total area.The second mutation occurs in the time range of 16.3-17.1µs, corresponding to the failure cloud map Figure 13c,d.Detonation products reach the end cover of the shell earliest at this stage.The crack is generated and expands inward along the outermost ring on the end cover.During this stage, the evolution time represents a mere 0.8% of the total time, while the growth area attributed to tensile damage accounts for 32.77% of the total area and the growth area associated with tensile fracture amounts to 18.71% of the total area.The two mutations account for 1.8% of the total time, while the damage and fracture areas account for 51.47% and 43.62% of the total area, respectively, indicating that the tensile failure of the shell mainly occurs in the middle and ends.
shear and tensile failure, respectively.The left side shows the cloud map depicting the failure face piece, while the right side shows the cloud map illustrating the presence of surface cracks on the shell.As can be seen from the cloud map, the shell crack propagates along the axial direction to both ends.Detonation products reach the end cover of the shell at 16.3 μs, completing more than 99.5% of the damage area growth.At 32.4 μs, the shear fracture area of 92.94% and tensile fracture area of 98.41% are increased, with a shell expansion for the initial diameter of 1.8 times.To analyze the discrepancies between naturally created fragments due to shear and tensile failure, the fractal of fragments at 60 µs is shown in Figure 14.At this specific duration the fracture of the shell is completed, leading to the enhanced visibility of fragments as they disperse uniformly.The detailed local magnification reveals three distinct regions.The two regions of the side wall centre and the end cover, correlating to the two mutations of tensile failure, exhibit lesser fragments volumes and fewer bar-shaped fragments relative to other areas.This suggests that tensile damage significantly influences the fragmentation outcome in both regions.As seen in the selected part of the dashed line frame in Figure 12b, the growth of the tensile damage/fracture area undergoes two mutations, with the mutation occurrence intervals for damage and fracture overlapping.The first mutation occurs in the time range of 5.5-6.5 μs, corresponding to the failure cloud map in Figure 13a.At this stage, the detonation products swiftly reach the shell wall, initiating the formation of a crack along the circumference starting from the center of the shell.During this stage, the evolution time represents only 1% of the total time, while the growth area attributed to tensile damage accounts for 18.70% of the total area and the growth area associated with tensile fracture amounts to 24.91% of the total area.The second mutation occurs in the time range of 16.3-17.1μs, corresponding to the failure cloud map Figure 13c,d.Detonation products reach the end cover of the shell earliest at this stage.The crack is generated and expands inward along the outermost ring on the end cover.During this stage, the evolution time represents a mere 0.8% of the total time, while the growth area attributed to tensile damage accounts for 32.77% of the total area and the growth area associated with tensile fracture amounts to 18.71% of the total area.The two mutations account for 1.8% of the total time, while the damage and fracture areas account for 51.47% and 43.62% of the total area, respectively, indicating that the tensile failure of the shell mainly occurs in the middle and ends.
To analyze the discrepancies between naturally created fragments due to shear and tensile failure, the fractal of fragments at 60 μs is shown in Figure 14.At this specific duration the fracture of the shell is completed, leading to the enhanced visibility of fragments as they disperse uniformly.The detailed local magnification reveals three distinct regions.The two regions of the side wall centre and the end cover, correlating to the two mutations of tensile failure, exhibit lesser fragments volumes and fewer bar-shaped fragments relative to other areas.This suggests that tensile damage significantly influences the fragmentation outcome in both regions.

Scattering Characteristics of the Natural Fragments
The displacement cloud diagram resulting from natural fragments at various moments is shown in Figure 15.This diagram represents fragments from both the side wall

Scattering Characteristics of the Natural Fragments
The displacement cloud diagram resulting from natural fragments at various moments is shown in Figure 15.This diagram represents fragments from both the side wall and end perspectives, as identified on the left and right sides, respectively.Notably, the displacement recorded at 100 µs is chiefly within the range of 12.62-16.84cm, by which time the acceleration of the fragments has stopped.To investigate the displacement characteristics of fragments and expansion characteristics of the shell at different acceleration stages, the dimensionless parameter β = Mag-Dis/RO is introduced to represent the relationships between the resulting displacement of the fragments and the shell radius.When calculating β at a certain moment, the average displacement value of all fragments at that moment is taken as MagDis.The relative spe- To investigate the displacement characteristics of fragments and expansion characteristics of the shell at different acceleration stages, the dimensionless parameter β = MagDis/R O is introduced to represent the relationships between the resulting displacement of the fragments and the shell radius.When calculating β at a certain moment, the average displacement value of all fragments at that moment is taken as MagDis.The relative specific volume can be considered as v = (β + 1) 3 .At v = 7 and v = 13; the value of β is 0.913 and 1.351, corresponding to 36.412 µs, 47.852 µs, respectively.
The resulting displacement distribution curves of fragments along the shell busbar are shown in Figure 16.These curves exhibit a smooth trend during the early scattering phase at 10 µs and 40 µs, transitioning to substantial fluctuations during the later scattering phase at 70 µs and 100 µs.This suggests that the dispersion of fragment displacement increases as the shell fractures and the fragments separate.At the corresponding moments of v = 7 and v = 13, the average resultant velocity of fragments is 3.844 cm and 5.688 cm, respectively.Utilizing the methodology presented in this paper improves the representation ability of fragment displacement by 47.97% compared to the cylinder test.The resulting displacement distribution curves of fragments along the shell busbar are shown in Figure 16.These curves exhibit a smooth trend during the early scattering phase at 10 μs and 40 μs, transitioning to substantial fluctuations during the later scattering phase at 70 μs and 100 μs.This suggests that the dispersion of fragment displacement increases as the shell fractures and the fragments separate.At the corresponding moments of v = 7 and v = 13, the average resultant velocity of fragments is 3.844 cm and 5.688 cm, respectively.Utilizing the methodology presented in this paper improves the representation ability of fragment displacement by 47.97% compared to the cylinder test.The fragments generated by shell fracture are accelerated by detonation products.Through the five marked points in Figure 8, the resulting velocity and displacement evolution of the location are recorded.Considering the symmetry between points P1 and P5 as well as that between P2 and P4, the evolution of the respective fragment parameters shares a convergent pattern.The fragments at P1, P2, and P3 are correspondingly labeled as Point1, Point2, and Point3.Figure 17 showcases the acceleration time history curves for these three fragments.The fragments generated by shell fracture are accelerated by detonation products.Through the five marked points in Figure 8, the resulting velocity and displacement evolution of the location are recorded.Considering the symmetry between points P1 and P5 as well as that between P2 and P4, the evolution of the respective fragment parameters shares a convergent pattern.The fragments at P1, P2, and P3 are correspondingly labeled as Point1, Point2, and Point3.Figure 17 showcases the acceleration time history curves for these three fragments.
Figure 17a shows the full curve of fragment acceleration.The fragments Point1, Point2, and Point3 are accelerated, with a final resultant velocity of 1.916 km/s, 2.034 km/s, and 2.086 km/s, respectively.Additionally, they attain 80% of their final velocities at 21 µs, 26 µs, and 32 µs.At the corresponding moments of v = 7 and v = 13, the average velocity of the fragments is 1.754 km/s and 1.847 km/s, respectively.The two values are 86.88% and 91.46% of the velocity at 100 µs, respectively.Utilizing the methodology presented in this paper improves the representation ability of fragment velocity by 5.30% compared to the cylinder test.Based on the analysis of displacement and velocity, the JWL equation of state can describe the scattering process of fragments over a wider range by calibrating it with DLCHEQ.
Figure 17b shows the initial curve of fragment acceleration.The earliest moment of fragment acceleration is at 4 µs and the latest is at 15 µs, which corresponds to the beginning of the fracture of the shell and the moment when the detonation products reach the end cover.This result shows that fragments develop from the centre of the busbar to both ends.In order to explore the acceleration law of global fragment scattering of the shell, The fragments generated by shell fracture are accelerated by detonation products.Through the five marked points in Figure 8, the resulting velocity and displacement evolution of the location are recorded.Considering the symmetry between points P1 and P5 as well as that between P2 and P4, the evolution of the respective fragment parameters shares a convergent pattern.The fragments at P1, P2, and P3 are correspondingly labeled as Point1, Point2, and Point3.Figure 17 showcases the acceleration time history curves for these three fragments.Figure 17a shows the full curve of fragment acceleration.The fragments Point1, Point2, and Point3 are accelerated, with a final resultant velocity of 1.916 km/s, 2.034 km/s, and 2.086 km/s, respectively.Additionally, they attain 80% of their final velocities at 21 μs, 26 μs, and 32 μs.At the corresponding moments of v = 7 and v = 13, the average velocity of the fragments is 1.754 km/s and 1.847 km/s, respectively.The two values are 86.88% and 91.46% of the velocity at 100 μs, respectively.Utilizing the methodology presented in this paper improves the representation ability of fragment velocity by 5.30% compared to the cylinder test.Based on the analysis of displacement and velocity, the JWL equation of state can describe the scattering process of fragments over a wider range by calibrating it with DLCHEQ.
Figure 17b shows the initial curve of fragment acceleration.The earliest moment of fragment acceleration is at 4 μs and the latest is at 15 μs, which corresponds to the beginning of the fracture of the shell and the moment when the detonation products reach the end cover.This result shows that fragments develop from the centre of the busbar to both ends.In order to explore the acceleration law of global fragment scattering of the shell, the resulting fragment velocity distribution curves along the shell busbar at 5 μs, 10 μs, 15 μs, 30 μs, and 100 μs are plotted in Figure 18.The velocity distribution of the resulting fragments can be divided into three stages.During stage A (0-18 μs) the shell undergoes partial fragmentation, experiencing acceleration from the center to both ends under the expansive force exerted by the detonation products.The resultant velocity of fragments at different positions is different.For example, at 15 μs the maximum velocity of the fragment reaches 1.529 km/s, while the end of fragment has a velocity of 0. During stage B (15-30 μs), detonation products are continuously loaded for the second time after the shell has completely broken.At 30 μs the maximum velocity of the fragment is 1.787 km/s, while its minimum velocity is 1.447 km/s.The average velocity of all fragments is 1.665 km/s, with 0.34 km/s separating the highest and lowest values, indicating that the velocity of fragments at different positions is close.During stage C (30-100 μs) there is a substantial dissipation of detonation products, leading to a decrease in the propelling force exerted by the fragments.At 100 μs the average resultant velocity of all fragments is 2.019 km/s, i.e., the average velocity of the fragments increases by 0.354 km/s in the time interval of 70 μs.
Stages A and B collectively constitute 30% of the total acceleration time, with the average velocity at 30 μs reaching 88.51% of the final average velocity achieved upon completion of the acceleration process.These stages serve as the primary acceleration stages for the fragments.The velocity distribution of the resulting fragments can be divided into three stages.During stage A (0-18 µs) the shell undergoes partial fragmentation, experiencing acceleration from the center to both ends under the expansive force exerted by the detonation products.The resultant velocity of fragments at different positions is different.For example, at 15 µs the maximum velocity of the fragment reaches 1.529 km/s, while the end of fragment has a velocity of 0. During stage B (15-30 µs), detonation products are continuously loaded for the second time after the shell has completely broken.At 30 µs the maximum velocity of the fragment is 1.787 km/s, while its minimum velocity is 1.447 km/s.The average velocity of all fragments is 1.665 km/s, with 0.34 km/s separating the highest and lowest values, indicating that the velocity of fragments at different positions is close.During stage C (30-100 µs) there is a substantial dissipation of detonation products, leading to a decrease in the propelling force exerted by the fragments.At 100 µs the average resultant velocity of all fragments is 2.019 km/s, i.e., the average velocity of the fragments increases by 0.354 km/s in the time interval of 70 µs.
Stages A and B collectively constitute 30% of the total acceleration time, with the average velocity at 30 µs reaching 88.51% of the final average velocity achieved upon completion of the acceleration process.These stages serve as the primary acceleration stages for the fragments.

Conclusions
This study analyzes the fracture and fragmentation of a thin-walled cylindrical shell under implosion force, and presents a strategy to gauge an explosive's work capability.The conclusions obtained are as follows: 1.
Using the thermochemical equation of state code DLCHEQ, the isentropic expansion curve of detonation products in the range of <13 was computed.Subsequently, the JWL EOS coefficients were calibrated across a wider range of temperature and pressure through curve fitting.This approach can successfully delineate an explosive's work capacity, as the derived equation of state accurately reflects the cylinder test's measurement results.Through this method, the fragmentation process of the shell can be described over a wider range.2.
Simulation of dynamic fracture process of the shell was undertaken based on the continuum-discontinuum element method.The validity of the numerical simulation results was corroborated through the Mott mass-quantity distribution theory and shelled charge explosion test velocity measurement.The discrepancy between the average initial velocity of fragments and the experimental value was 9.04%.The mass-quantity distribution curve of the fragments aligned well with the theoretical curve, demonstrating a discrepancy of only 5.83% in the total fragment count when compared to Mott's theoretical prediction.From these results, it can be seen that the scattering process of natural fragments generated by the thin-walled shell under implosion load is accurately described by the method.

3.
The thin-walled shell employed in this study primarily failed through shearing, generating strip fragments.In the evolution of tensile failure, two mutations generate local particle fragments, which constitute 1.8% of the total time and cover over 50% of the total failure area.Compared to shear failure, there is a "hysteresis effect" at the beginning of tensile failure and the hysteresis period is 0.5 µs.The shell's failure process mainly lasts for 4-18 µs, after which the unit's damage area stops growing and the fracture area begins to expand slowly.According to the processes of shell fracture and fragmentation, the acceleration process of natural fragments can be divided into three stages.The primary stages of fragment generation and acceleration are stages A and B, respectively.At 30 µs the shear/tensile damage/fracture area accounts for over 90% of the final value, and the average resultant velocity of fragments is 88.51% of the final value.The velocity and shear fracture area continue to gradually rise in stage C, indicating that fragments continue to be accelerated by detonation products outside the calibration range of the cylinder test.The final resultant velocity of fragments ranging from 1.8 to 2.1 km/s at 100 µs.

Figure 1 .
Figure 1.Schematic diagram of EOS of detonation products using DLCHEQ.

Figure 4 .
Figure 4. Configuration of cylinder test in explosive tower: (a) global view and (b) physical view.

Figure 5 .
Figure 5. Results of the cylinder wall expansion curve test and calculation.The accuracy of the JWL EOS was verified over the range of relatively specific volumes measured in the cylinder test.The 25 mm cylinder test can describe the work capac-

Figure 4 .
Figure 4. Configuration of cylinder test in explosive tower: (a) global view and (b) physical view.

Figure 4 .
Figure 4. Configuration of cylinder test in explosive tower: (a) global view and (b) physical view.

Figure 5 .
Figure 5. Results of the cylinder wall expansion curve test and calculation.The accuracy of the JWL EOS was verified over the range of relatively specific volumes measured in the cylinder test.The 25 mm cylinder test can describe the work capac-

Figure 5 .
Figure 5. Results of the cylinder wall expansion curve test and calculation.
Numerical model (b) Block model (c) Interface model

Figure 8 .
Figure 8. Schematic diagram of the validating model of the shelled charge.

Figure 8 .
Figure 8. Schematic diagram of the validating model of the shelled charge.

Figure 9 .
Figure 9. Schematic diagram of structural surface softening model: (a) tensile failure and (b) shear failure.

Figure 9 .
Figure 9. Schematic diagram of structural surface softening model: (a) tensile failure and (b) shear failure.

Materials 2023 ,
16,  x FOR PEER REVIEW 16 of 24 of milliseconds.The fracture and fragmentation process of the shell were analyzed based on the verified numerical simulation results in the CDEM program.

Figure 12 .
Figure 12.Damage/fracture spring area of the element under tensile/shear forces: (a) shear damage/fracture area of spring and (b) tensile damage/fracture area of spring.

Figure 12 .
Figure 12.Damage/fracture spring area of the element under tensile/shear forces: (a) shear damage/fracture area of spring and (b) tensile damage/fracture area of spring.

Figure 17 .
Figure 17.Velocity-time history curve of fragments in acceleration stage: (a) full acceleration curve and (b) starting moment of acceleration.

Figure 17 .
Figure 17.Velocity-time history curve of fragments in acceleration stage: (a) full acceleration curve and (b) starting moment of acceleration.Figure 17.Velocity-time history curve of fragments in acceleration stage: (a) full acceleration curve and (b) starting moment of acceleration.

Figure 17 .
Figure 17.Velocity-time history curve of fragments in acceleration stage: (a) full acceleration curve and (b) starting moment of acceleration.Figure 17.Velocity-time history curve of fragments in acceleration stage: (a) full acceleration curve and (b) starting moment of acceleration.

Table 1 .
Coefficients of the JWL EOS from calibration.

Table 1 .
Coefficients of the JWL EOS from calibration.

Table 2 .
Velocimetric target data at 4 m and 5.2 m.

Table 3 .
Coefficients of the JWL EOS from calibration.

Table 4 .
: Distribution of fragments number at different quality intervals.

Table 4 .
Distribution of fragments number at different quality intervals.