Impact of Imperfect Kolsky Bar Experiments Across Different Scales Assessed Using Finite Elements

Kolsky bars or split Hopkinson pressure bars (SHPB) have been growing in popularity for testing materials at high rates since the 1990s [1–5]. A SHPB system consists of 3 long rods usually made of the same material and axially aligned called the striker (SB), incident (IB), and transmitted (TB) bars. In a compression test, a right cylindrical sample is placed between the IB and TB with the contact surfaces well lubricated, at which point the SB is launched with compressed air to impact the free end of the IB. This interaction causes a stress wave to propagate down the IB and interact with the sample, where a portion of the wave is reflected back into the IB and a portion moves through the sample and transferred to the TB. A strain gauge placed at the center of the IB and TB is used to measure these waves from which the response of the sample can be determined. The process and system are described in great detail in the Refs. [1–28].

Most systems are in the 10–25 mm diameter range and test at rates of 10²–10³ s⁻¹ [1–17]. Smaller bars with the accompanying smaller diameters are needed to push samples to higher strain rates of the order of 10⁴ and 10⁵ strain per second [18–21,23]. Some finite element (FE) modeling of Kolsky systems has been conducted, but the detail involved varies from publication to publication. The fidelity of the modeling techniques and the level of detail presented in these models can make it difficult to reproduce these simulation [19,24–27]. Simulations typically assume perfect alignment and do not model bar support points or gravitational effects [19,24–27]. By running several simulations modeling various misalignments or imperfections, the uncertainty in the output of such experiments can be quantified directly. Such an investigation has not been presented in the literature to date, so the uncertainty quantification conducted with the use of these simulations can aid researchers in understanding the variability of experiments and impact of these potential sources of error.

We provide an easy to follow modeling approach, describing in detail how to properly tune a simulation for any length and diameter scale. We demonstrate the usefulness of highly optimized symmetry-based simulations for mesh convergence and temporal resolution studies, highlighting the need for higher sampling frequencies to accurately capture the measured pulses in smaller diameter systems compared to larger ones. This is followed by the creation of a three-dimensional, half symmetry model, which is used to simulate real-world defects or imperfections. Specifically, we examine the effects of gravity and friction on the contact between the bars and their bushing supports and the effect of misalignment of the bars on the predicted material response of an aluminum 2024 sample. Simulations are conducted on a 12.7-mm diameter and 3.16-mm diameter system to evaluate the potential change in the impact of these effects across two different length scales.

Before we describe our simulation process, we first define three different materials which will be used in the simulation. First is C350 Steel which will be used when modeling the striker, incident, and transmission bars. This material has an elastic modulus of 202 GPa, a density of 8091 kg/m³, a Poisson's ratio of 0.272, and a yield strength of approximately 2 GPa. Brass is used for the bushings or collars and has a 8400 kg/m³ density, an elastic modulus of 117 GPa, and a Poisson's ratio of 0.34. Finally, the sample is modeled as 2024 aluminum using the Johnson-Cook (JC) plasticity model presented by Millan et al. and shown in Table 1 [28].

The first set of simulations were axis-symmetric and were good tools for determining features such as mesh density and temporal sampling frequency in explicit dynamic simulations. While mesh density is usually a well-understood concept, temporal frequency is something unique to explicit simulations. These Kolsky experiments occur on time scales of hundreds of microseconds and the data capture rate can greatly impact the solution, especially in unpulse shaped Kolsky bar simulations in which high-frequency content is present. A 3.16 mm diameter and 12.7 mm diameter system are modeled to capture the miniature and standard sizes of experiments and to reveal any differences that misalignment or gravitational effects may have on simulation results. The 3.16 mm system has a striker of 99.16 mm in length and incident and transmission bars 812.8 mm long, matching the miniature system used at Penn State. The 12.7 mm system has a 146.4 mm long striker and main bars with a length of 1.2 m. The modeling process will now be described, and is visually represented in Fig. 1.

We begin the process with the smaller diameter Kolsky system, as the small bar is expected to contain a higher frequency signal. This will also necessitate a smaller element size along the length of the bars. To determine the appropriate temporal frequency, we start with an overconservative mesh of 0.2 mm long elements with 8 elements across the radius and run simulations sampling at 2.5, 2.0, 1.5, and 1 MHz to determine the appropriate converged sample rate. This results in a 2 MHz sampling rate for this small bar system as shown in Fig. 2, which depicts the centroidal output of the element on the surface of the bar and closest to the middle of the incident bar, which represents the location of a strain gauge in a physical system.

With the proper temporal frequency identified, we filter the results moving forward using a lowpass filter and a cutoff frequency of ¼ of the sampling rate. This will prevent aliasing effects as only a signal of roughly ¼ the sampling rate can be accurately captured.

Next, we decrease the mesh density in stages. We begin with increasing the length of the elements to find the converged length. When identified, we then reduce the number of elements in the through thickness. Finally, we run the same procedure on only the striker to further reduce its mesh density. Since the striker's only role is to input the stress wave to the rest of the system it does not need to carry a highly accurate wave itself, but rather transfer a wave to the incident bar accurately. The convergence plots for each step as well as the final vs. the original mesh are shown in Fig. 3 with the selected plots in black.

Once the mesh density for the bars was set, more simulations were conducted using an aluminum 2024 sample with properties provided in Sec. 2.1. The sample geometry is based on a length/diameter ratio of 0.8 and a sample diameter to bar diameter ratio of also 0.8. Both the incident and transmitted waves are analyzed, just as in a real-world experiment, and these pulses are used to determine when the sample mesh is properly refined. The results of this refinement process are shown in Fig. 4, with the selected mesh results again in black. The coarsest mesh used, 8 × 0.12 mm has nearly identical performance without degradation of the wave behavior.

For the contact modeling, two different property sets were used. Set one had a hard normal behavior and a frictional coefficient of 0.1 while set two had a hard normal contact behavior and a frictional coefficient of 0.05. Set one was used for the bar on bar contact of the striker on the incident bar and set two was used for the bar/sample interfaces. This low frictional coefficient replicates a well lubricated surface, which is a standard practice in the experimental field. The striker/incident bar surface to surface interaction used the penalty contact formulation. This formulation was not stable for the direct contact between the sample and bars, so the kinematic contact formulation was implemented.

The same process of temporal and mesh convergence was followed for modeling the 12.7 mm diameter system. The final parameters for each model are given in Table 2, and incident, reflected, and transmitted pulses for the converged 12.7 mm system are shown in Fig. 5.

For our modeling of misalignment and gravity, we used a half symmetry model with the plane of symmetry being x − y. These models were meshed using the mesh density parameters and frequency of sampling found with the axis-symmetric simulations. Two times the number of elements along the radius are meshed along the half circumference. For example, 12 elements are used for the half symmetry of the small bars. Additionally, the parts are meshed with a hex-dominated structure such that the center core of elements are tetrahedral elements instead of hex elements. This was done to maintain symmetry in the mesh, which is reflective of the geometry. An example of the meshing for the 12.7 mm case is shown in Fig. 6, which shows a zoomed-in section centered on the sample. The 3.16-mm case has nearly identical element quality and is visually very similar to the that shown in Fig. 6.

Bushings or collars are also added to this model, replicating supports for the main bars, which are present in experiments. The bushing sleeves are modeled using brass material properties presented in Sec. 2.1. The outer radius is 11.4 mm, and the inner radius is 50 micrometers larger than the respective bar radius. The bushing was meshed coarser than the underlying bar by approximately 50% along the length, two fewer elements along the circumferential direction, and three and two elements along the thickness for the small and large systems, respectively. The bushings act as little more than a source of support and friction for the bars, so their meshing does not require high levels of accuracy. Additionally, having the higher density brass section also meshed coarser than the bar allows it to be assigned as a pure master surface in the contact algorithm. These interactions were modeled using interaction set two described in Sec. 2.2 using the kinematic contact formulation. Four bushings are placed on each bar. The center of the first and last bushings is 2 cm from the ends of the bar, with the other two equally spaced in the middle.

To model the imperfection and the dynamic impact event accurately, we broke the simulation into pieces, a static general step to solve for the deformed configuration caused by imperfection and an explicit dynamic step for the impact event and resulting wave dynamics. To accomplish this, we requested a restart file be generated at the end of the static step, which is used to create a predefined field in the dynamic simulation. As a result, the separate static simulations need to be run to completion first. This differs from a standard multistep process because of the change from the implicit static general solution step to the explicit dynamic solution step. For all models, the bars began perfectly aligned with the striker and sample at the beginning of the static step. A flowchart for the modeling process is shown in Fig. 7.

For the gravity model, the incident and transmission bars start in contact with the bushings. The outsides of the bushings were locked in all directions and a gravity load was placed on these bars only. The striker was not included because it has no such supports in an experiment and the sample was not included because the low friction surface prevented a converged solution. The mass of the sample is so low relative to the rest of the system that any effect from the sample-related gravity forces would be negligible.

The misalignment was simulated in four different misalignment cases by displacing certain bushings in the static step. The bushings were still locked in rotation as well as x and z displacements, but a y displacement condition was prescribed and the bushings began in contact with the bar at the appropriate side. For each case, the bushing closest to the striker was not moved and the last bushing on the incident bar and first bushing on the transmitted bar was moved by the same amount in the same direction. This was done because the striker/incident and incident/transmission bar interfaces can be very precisely aligned with each other experimentally, so the relative alignment in simulations is assumed to be perfect. The four cases are described by the bushing displacements in Table 3. These cases of alternating displacements positively and negatively and cases of displacing the edges down with the center up create the worst-case bending conditions in the bars given a fixed displacement magnitude. By taking data from elements on the top and side of the bar, still in the center length-wise, we can collect information for a horizontal and vertical misalignment case from the same simulation.

It should be noted that due to the finer mesh, the smaller bar took longer to run the static and dynamic steps and the static cases did have trouble converging. The initial step size for the misalignment cases had to be adjusted down to allow the solution algorithm to proceed. Once the first 1–2% of the static simulation was completed, the process sped up quickly with the step size increase for each following increment of load. This is the result of the contact algorithm struggling to equilibrate and generate the deformed shape of the bars. Once that was done, applying further load was trivial as evidenced by a steady increase in step size for the remainder of the static simulation.

Each static simulation was run using 10 cores and 2 GB of RAM each, with each dynamic simulation running on 20 cores with 2 GB of RAM each. The static simulation run times varied from 15 min for some of the larger bar diameters simulations to about an hour and a half for the longest small bar simulation. The large bar dynamic simulations took about 45 min to run while the smaller bar dynamic simulations took about 1.6 h. This is significantly longer than the 1–3 min it took the axis-symmetric simulations to run on a single core, which shows the effectiveness of such small simulations for large throughput, quick turnaround simulations for activities such as mesh refinement. That being said, the axis-symmetric model results did fall short compared to an ideal half symmetry simulation as can be seen in Fig. 8. The IB pulses, not shown, have the same slope and peak, with later time values being overpredicted by the axis-symmetric model by 1–2%. As shown in Fig. 8, the reflected and transmitted pulses for both length scales do not have the proper magnitude in the axis-symmetric model but do have the right profiles. This is believed to be due to the lack of accurate radial effects captured in the axis-symmetric modeling case compared to the half symmetry case using three-dimensional elements. The difference presenting itself in the reflected and transmitted pulses also indicates that the sample interactions and plastic deformation are not accurately reported in an axis-symmetric simulation.

Figures 9 and 10 show the incident, reflected, and transmitted pulses for each imperfection accompanied by the ideal half symmetry model for reference.

Keep in mind that while the incident pulse is a measured quantity and is used to evaluate system response and performance, it acts only to load the sample, and therefore it is not as important if the incident pulse changes as a result of these imperfections. The reflected and transmitted pulses are directly related to a sample's strain rate and stress during the test, making the analysis of these waves more important.

For both size scales, misalignment cases 1&2 and 3&4 overlap for the incident and reflected pulses while the transmitted pulse appears largely unchanged between the misalignment cases. It should also be noted that horizontal misalignment did not produce significantly different results from one another, due to minimal stresses being generated in the locations of interest. The misalignment cases did differ from the gravity simulation for the large bar, but the gravity and reference configurations do not differ to the same degree as in the small bar case where the gravity and reference configurations are distinct.

In the small bar, the incident pulse is not greatly impacted by any of the imperfection cases with the only changes visible in later time. The application of gravity caused a pre-existing positive stress in the bar for the reflected wave, which can be seen in the 240–248 microsecond period of Fig. 9 of 5–15 MPa of pre-existing stress in the measurement, which is maintained for most of the time duration of the pulse. For the horizontal misalignment, all of the profiles fell in line with the reference case as can be seen in the middle right plot of Fig. 9. It can also be clearly seen that gravity causes the latter portion of the reflected pulse to level out while the remainder of the cases continue to decline with time. This difference equates to upwards of a 50% increase in the stress profile in the 275–283 microsecond time frame compared to the reference case. There is a significant spread in the vertical misalignment cases in the 260–280 microsecond period, showing the sensitivity in magnitude and profile of the smaller system to bushing misalignment. This difference is as high as 36% variation from the average depending on the misalignment case.

The large bar data are less black and white than those for the smaller bar. The horizontal misalignments continue to cluster together despite containing different load cases, but they are all distinctly larger in magnitude for the reflected pulse than the reference or gravity cases, constituting a 10–20 MPa increase over the majority of the pulse duration from 370 to 420 microseconds. For vertical misalignments, the 1&2 load cases are distinctly separate from the 3&4 and gravity loads in the incident pulse in the 150–180 microsecond range. The vertical misalignments in the reflected pulse are fairly identical until the unloading section of the wave around 420 microseconds, indicating a lesser sensitivity to misalignment errors. There is roughly a 5–20 MPa variation from the reference configuration. This variation is almost exclusively positive relative to the ideal reference case, which differs from the small bar case in which the misalignments produced both positive and negative variations. The large bar has a difference of approximately 18% variation from the average across the reflected pulse, which is half the variation seen in the smaller bar. We also see a minimal difference between the reference case and the gravity simulation, especially when compared to the small bar system. This is likely due to the higher input of energy to the 12.7-mm diameter system from the larger and longer striker compared to the striker in the 3.16-mm system. This larger relative amount of energy means losses due to friction have a lower overall impact on the wave dynamics. This result means that small bars are also more sensitive to gravitation loads.

After a careful and intensive model generation process, we have produced a rigorous half symmetry model using Abaqus® for two different size scales of Kolsky bar systems. We demonstrated the higher fidelity of these models compared to simpler axis-symmetric simulations, and show how our half symmetry system can be used to model imperfections that will inevitably be present in real-world experiments. We clearly show the sensitivity of smaller systems to bushing misalignment, which manifests itself in the reflected pulse and therefore directly impacts the strain rate seen by the sample. Even though the larger system has more mass and therefore the potential to be more impacted by the simulation of gravity, we see that the smaller bar is more affected. We then ultimately conclude that smaller Kolsky bars are more sensitive to imperfections overall compared to their larger cousins.

A limitation of the existing work is the lack of ability to precisely offset the bushings in the 50-μm increments defined. This value was selected to replicate real-world tolerancing issues with such systems, making it not feasible to conduct such experiments with reproducible accuracy at this time. Future work will entail gathering experimental data, as well as procuring hardware and modifying existing experimental equipment to be able to prescribe such minute offsets. Additionally, several experiments of a perfectly aligned system will be conducted to establish a statistically significant response to which this simulation structure can be compared.

Los Alamos National Laboratory (Subcontract No. 424285; Funder ID: 10.13039/100008902).

The authors attest that all data for this study are included in the paper.