ABSTRACT
Stereoscopic digital image correlation (DIC) is used to measure the shape evolution of a soft, transparent thermoplastic elastomer subject to a high strain rate compression test performed using a Kolsky bar. Rather than using the usual Kolsky bar wave analysis methods to determine the specimen response, however, the response is instead determined by an inverse method. The test is modeled using finite elements, and the elastomer stiffness giving the best match with the shape and force history data is identified by performing iterative simulations. The advantage of this approach is that force equilibrium in the specimen is not required, and friction effects, which are difficult to eliminate experimentally, can be accounted for. The thermoplastic is modeled as a hyperelastic material, and the identified dynamic compressive (non-linear) stiffness is compared to its quasi-static compressive (non-linear) stiffness to determine rate sensitivity.
INTRODUCTION
Tissue simulant materials are sought to provide realistic experimental devices to simulate the human body’s response to blast or impact loading that can occur in military scenarios, law enforcement and emergency response events, vehicle accidents or sporting events [1]. This approach is meant to help develop better protective equipment or procedures to prevent serious injury or death. In most practical injury scenarios, tissues are subject to dynamic loading involving large amplitude strains (20 % [2]) to vulnerable soft tissues at strain rates above 10 s-1 [3]. Numerical models of these test devices are crucial to interpreting the measurement data in these complicated tests, and efforts are underway to provide the material data needed to calibrate such models. Simpler uniaxial mechanical tests of soft tissue specimens show that the large-strain response of these materials is generally non-linear and rubber-like, and can be represented by hyperelastic models developed for polymers [2, 4, 5]. Strain rate sensitivity is also generally observed in soft tissues, prompting the use of viscoelastic models to describe the relaxation behavior [2, 6, 7].
High strain rate measurements of actual tissues and tissue simulants are often performed using a Split Hopkinson Pressure Bar, or Kolsky Bar [8, 9], which can achieve large strains at uniform strain rates in excess of 100 s-1. The difficulties in obtaining valid high strain rate data on soft materials using Kolsky bar techniques are well documented. Soft materials take much longer to achieve mechanical equilibrium when subject to a rapidly changing load. In most practical situations, equilibrium is not established in the sample until the test is nearly over, invalidating most if not all of the test [8]. Careful selection of specimen thickness and the use of pulse shaping to increase the rise time of the load pulse have been effective methods for achieving valid dynamic test results [10]. Measuring the forces on soft samples is also challenging because they are below the typical sensitivity of Kolsky bars designed for testing metals. Special, highly sensitive force transducers placed directly on either side of the sample are required to obtain force signals [11]. Care must be taken to separate out inertial effects in these force signals to obtain the true specimen response [12]. Finally, confinement techniques have been successfully employed to force either hydrostatic or shear loading conditions at high strain rates [13].
For extremely soft materials (less than 10 Shore A), such as the material examined in our experiments, achieving force equilibrium in a Kolsky bar test is exceedingly difficult. Further, friction between the sample and the compression platens is difficult to eliminate entirely, even with generous lubrication. Under dynamic loading, these effects will tend to produce non-uniform deformation in the sample in the form of large amplitude surface waves and prominent barreling which invalidates the usual Kolsky bar assumptions of uniform, uniaxial stress and strain throughout the specimen. In low-rate compression tests, the effect of friction on the perceived stress-strain behavior can be significant. This has led some researchers to impose no-slip boundary conditions in order to more accurately determine the material response with a known, rather than ambiguous, state of friction [4, 14].
Recent advances in optical shape measurement using stereoscopic (3D) Digital Image Correlation (DIC) [15] and in high speed digital camera technology has now made it feasible to measure the shape evolution of soft specimens during dynamic testing in a Kolsky bar. Further, by combining this new measurement capability with finite element modeling, one may be able to deduce the constitutive behavior of the material at high rates using so-called inverse methods [16]. Such methods, which include minimizing the error of finite element simulations against experimental data by systematically adjusting the relevant material parameters has been used to identify the properties of metals [17-19], composites [20], ceramics [21], polymers [22] and biomaterials [23] usually under quasi-static test conditions but more recently at high rates of strain [24].
In this paper, high-speed 3D DIC is used to measure the dynamic deformation of a soft tissue simulant material while it is compressed at high strain rate using a modified Kolsky bar technique. The sample is allowed to deform non-uniformly due to specimen ring-up effects and to friction between the sample and the bar. Finite element modeling is employed to deduce the constitutive behavior of the material by systematically varying the physical parameters governing the mechanical behavior of the modeled specimen and the friction coefficient to match the DIC history and force history data. This approach was first described in [25]. In this updated work, commercial software is used to perform model sensitivity analysis and to identify optimal parameter values. The identified optimal parameters are averaged for five individual experiments and the result is then compared to the quasi-static behavior to determine rate sensitivity of the material. Standard deviations of the identified parameters are used to estimate the overall uncertainty of the identified dynamic response of the specimen.
EXPERIMENTAL METHODS
The elastomer was cast into 9.5 mm diameter by 6.5 mm thick specimens for dynamic compression testing. The material has a density of 870 kg/m3, and it is nearly transparent. Tests were conducted using a maraging steel Kolsky bar measuring 15 mm in diameter with bar lengths of 1500 mm. Forces were measured with piezo-electric dynamic force transducers placed on either side of the specimen, as shown in Fig. 1 below. The force transducers have a resolution of ± 0.1 % of full scale, or ± 0.5 N on the incident side and ± 0.25 N on the transmission side. Polished steel platens (5 mm thick by 15 mm diameter) are placed between the transducers and the sample. No lubrication was used in the experiments. Compressive pulses were created by impacting a 250 mm long, 15 mm diameter maraging steel striker bar against the incident bar at 5.5 m/s.
An additional difference between these tests and more traditional Kolsky bar testing is that the upstream platen and force transducer are allowed to collectively decouple from the incident bar to become a "flyer plate." The flyer has the advantage of producing a much larger dynamic strain in the specimen than if the transducer and platen remained fixed to the incident bar, which was desirable given the limited compression pulse width (and thus specimen strain) our Kolsky bar could generate on a sample of this size. The flyer test is accomplished by using a low-tack adhesive at the incident bar-transducer interface that is strong enough to hold the pieces in place prior to the test but too weak to resist the tensile wave reflection off the incident platen as the wave passes through the sample. A thick, soft rubber pulse shaper was employed to produce a long, gradually rising stress pulse, which helped develop a repeatable release of the flyer from the incident bar.
DIC measurements of the sample deformation are acquired at 120,000 frames per second with an image size of 128 x 208 pixels. Commercial DIC software1 [26] is used to perform the image correlation measurements. The cameras use 90 mm macro (1:1 magnification) lenses with f/5.6 apertures and are placed 30 cm from the sample with a 12.5 ° pan angle. The resolution is 18 pixels/mm.
Fig.1 Axisymmetric finite element model showing the arrangement of the sample, platens, load cells and transmission bar
A random speckle pattern was created on the samples using a light dusting of flat black spray paint. Attempts were made to achieve a nominal average speckle size of 5 to 7 pixels in diameter and a coverage factor of 50 % following manufacturer’s recommendations. The speckle patterns were backlit using a Teflon reflector placed directly behind the specimen and illuminated by halogen optical fiber lamps, as shown in Fig. 2. Using this optical set-up, rigid body translations were measurable within ± 0.01 mm. Fig. 3 indicates the typical correlation region on a sample. Correlation measurements used the default software analysis settings (21 pixel windowing, 5 pixel overlap, with default smoothing). The DIC measurement data, obtained from about 200 stereo image pairs during the test, are automatically oriented to a reference plane that is defined by fitting the initial shape of the (cylindrical) specimen prior to deformation. Because of this automatic plane fit and the fact that the sample deforms axisymmetrically, the DIC coordinate system could be aligned with the axisymmetric finite element model coordinate system using simple y- and z-translations.
FINITE ELEMENT MODEL
The finite element model consists of the flyer assembly (force transducer plus steel platen), the sample, the transmitted force transducer and platen, and the transmission elastic bar. ABAQUS/Explicit1 is used to perform the simulations. A portion of the modeled mesh was shown earlier in Fig. 1. The model uses axisymmetric CAX4R elements, with 352 elements in the sample and much coarser meshes in the steel parts. The maraging steel compression bars as well as the platens are modeled as linear elastic solids with the following properties: E = 1.9×1011 Pa, v = 0.29, and p = 8048 kg/m3. The donut-shaped force transducers are modeled with the same stiffness and Poisson ratio as the maraging steel but with a density of 6594 kg/m3. The boundary condition for the simulation is set by specifying the velocity history of the flyer, which is obtained directly from DIC measurements on the flyer itself. Simulations are carried out until the flyer velocity vector begins to deviate from the compression axis, or "pitch," which occurs eventually when the flyer decelerates and comes to rest against the specimen.
The elastomer is modeled as an isotropic hyper-elastic material using the Marlow strain energy potential [27]. Hyper-elasticity is characterized by large, recoverable elastic strains that are characteristic of rubbery materials. Rubbery response in polymeric materials occurs in the region nestled between the visco-elastic response and viscous flow region in strain-rate space [28]. Micromechanical behavior associated with rate dependency and dissipation (viscoelastic effects), such as chain slippage and the breakage and reforming of secondary (crosslink) bonds, is relatively unimportant in the rubbery regime.
Fig. 2 Measurement set-up
Fig. 3 Typical correlation region shown on an actual image of a speckled sample from which correlation measurements are obtained
The rate sensitivity of the elastomer is examined by attempting to identify a purely hyperelastic dynamic response by inverse methods and comparing this identified response to the quasi-static behavior. If the identified dynamic response is found to be much different than the quasi-static response, rate sensitivity would thereby be indicated. The baseline constitutive response of the material used in the model is derived from uniaxial stress-strain data obtained at a low strain rate (4 x 10 -4 s-1) following ASTM D 575-91 [29]. The quasi-static response, shown in Fig. 4, is characteristic of rubbery materials, showing an upturn in the stress-strain curve due to the finite deformability of even very long chain polymeric molecules. To maintain this basic character while allowing for possible rate effects, the response of the specimen is modified by a stiffness scaling factor, m, as shown in Fig. 4.
Fig. 4 Quasi-static compression stress-strain response of tissue simulant with polynomial fit and two alternate models with twice (m = 2) and half (m = 0.5) the stiffness of the quasi-static response (defined as m = 1)
Rayliegh damping is employed in the numerical scheme to eliminate unphysical oscillations in the simulations to better mimic the dynamic response of the sample. This method is employed by the finite element code as a generic means to account for dissipation in many different materials [27]. It works by adding a small damping stress, ad, to the stress from the basic hyperelastic response that is proportional to the strain rate, by specifying a positive damping factor fi:
The magnitude of the damping factor must be chosen with caution as it is not intended to simulate the strain rate effects due to the bulk micromechanical response of the specimen. To avoid obscuring gross strain-rate effects from those related to numerical damping, a minimal damping coefficient must be chosen that is just large enough to eliminate non-physical oscillations while not adding significant numerical stiffness to the sample.
A commercial software package [30]1 is used to perform the model sensitivity analysis and to identify the dynamic sample stiffness by comparing the simulation results to the experimental data. The software acts as GUI-driven macro that alters finite element model input data, controls the solver execution and displays and analyses the simulation results. It also approximates the finite element model response over the variable space of interest using interpolation functions, and employs optimization tools to identify optimal parameter values using objective functions describing the agreement between the simulation results and the data. The objective functions are described next.
OBJECTIVE FUNCTIONS
Objective functions expressing the difference between the finite element model results and the data are built individually for the force history and shape history data. Then, a sing le "cost" function is assembled to represent the overall agreement between simulation and data for identifying optimum parameter values. DIC shape history data are compared to the simulation nodal displacements as indicated in Fig. 1. Residuals are computed at the DIC measurement locations along the center of the correlation region parallel to the compression (y-)axis at each measurement time point (e.g. for each image pair acquired and analyzed). The modeled surface positions are interpolated to match the DIC measurement positions along the length of the specimen. The shape history objective function, 0S, is given by:
In Equation 2, Z is the out-of-plane surface position of the deforming sample. The superscript FEA refers to the finite element model result, while the superscript DIC denotes experimental data (DIC measurement). Subscripts i and t refer to space and time, respectively. The transmitted force history objective function, 0F, is given by:
F is the transmission force, while EXP refers to the experimental data and FEA is again simulation results. At each time step, 0F is calculated by linearly interpolating the Ffea at the locations where FEXP is available. The combined objective function, 0, is assembled as follows:
In the above equation, SFf and SFS are scale factors for the force and shape residuals, respectively. The scale factors are selected to weight each individual objective function such that the order of the scaled objective is 1 near the optimum point. For this study, SFf = 1000 N and SFS = 0.001.
RESULTS
Force-time data from a typical experiment are shown in Fig. 5 along with a sketch showing the corresponding behavior of the flyer plate during the test. Prominent features of the data are labeled in the graph. A sharp rise in the incident force marks the arrival of the incident compressive strain wave generated by the striker. As this wave reflects from the incident platen-sample interface, a tensile (negative) force is observed. This tensile load causes the bond holding the flyer to fail, releasing it into the specimen. As the specimen compresses further, the transmitted force steadily increases until the flyer arrests, causing the transmitted force to peak. Next, the sample begins to release its stored strain energy, pushing back on the flyer. Soon the flyer is thrown back against the incident bar, which all the while has been advancing forward due to the action of the trapped compression wave. This second impact is marked by the final sudden rise in both force signals. Because the free surface of the flyer carries no stress, the forces measured on either side of the specimen will never be equal. Because of this, the force signals cannot be compared to check sample equilibrium nor can they be relied upon to determine the sample stress in the usual way.
Fig. 5 Force signals recorded on the incident and transmitted side of the specimen during a flyer experiment (left) and a sketch of the flyer experiment as visualized using the axisymmetric finite element simulation (right)
Fig. 6 plots the overall sample strain and strain-rate versus time for the same experiment. In the test shown, almost 70 % engineering strain is achieved before symmetry breaks down. A Kolsky bar 10 times as long as the one used here would be needed to achieve this level of strain using the same sample at this strain rate without using the flyer technique. A second observation is that the strain rate is relatively uniform over most of the test. Ordinarily this is critical in normal Kolsky bar tests. Here, however, it is less so in this study because the assumption of uniform strain rate is not used in the inverse analysis to deduce the specimen response.
The finite element model parameters governing the mechanical response of the sample are the stiffness factor, m, the friction coefficient, f, Poisson’s ratio, v, and the damping factor, p. In principle, all of these factors can be examined simultaneously by conducting a sensitivity analysis using a large Design-of-Experiments (DOE) matrix. However, certain parameters can be specified ahead of time to reduce the number of unknowns in the problem. For example, since very soft polymeric materials are incompressible, v = 0.5 should be prescribed. However, computational stability requires that some compressibility be added, which will affect accuracy of the solution of this highly confined compression test [27]. Computational cost and accuracy were found to be adequate with v = 0.4950. The penalty for allowing compressibility is a systematic under-prediction of the out-of-plane deformation, which introduces a systematic error in the value of 0S. Another parameter that must be prescribed is the damping factor, p. Selecting p = 0.000025 s provides realistic-looking force signal while having little overall effect on the force levels themselves, and therefore not adding a gross amount of stiffness to the material that would confuse attempts to identify a dynamic value of m. Using Eq. 1, this damping level adds a stress equal to about 1 % of the zero-strain modulus of the material for a global strain rate of 400 s-1.
Fig. 6 Global strain and strain rate histories for a typical flyer experiment from DIC-measured flyer velocity
With v and p fixed, we proceed to investigate the influence of the stiffness factor, m, and the friction coefficient, f. To this end, a Full-Factorial DOE is performed between 0.0 < f < 2.0 and 0.75 < m < 1.25 with 10 levels for each parameter. The influences of stiffness and friction on 0F and 0S are indicated in Fig. 7. The friction coefficient has dramatic influences on both 0F and 0S at low values, but at higher values the sensitivity to friction is minimal. For 0F, the sensitivity to friction falls beyond about 0.3, while for 0S the influence dramatically lessens above 0.5. Due to the lack of sensitivity above f = 0.5, identifying an optimal friction coefficient would be difficult and the result would likely be unreliable. Consequently, the friction coefficient is set to the limiting case of no-slip for the remaining simulations used to identify m.
DOEs were executed for each experiment with 0.75 < m < 1.5 to identify the dynamic stiffness of the specimen relative to its quasi-static response. The other parameters are fixed to values discussed previously, namely: p = 0.000025 s, v = 0.495 and no-slip friction. A representative plot of the effect of m on 0F and 0S is shown in Fig. 8. Both 0F and 0S have distinct minima, though not at the same value of m. A conclusion from the latter observation is that the model is not able to achieve perfect agreement with experimental observation, leading to this tradeoff between shape and force agreement. Since neither objective has a clear precedence over the other, the objectives will maintain approximately equal weighting by using the previously-defined scale factors. The simulation results of Fig. 8 are fit with a Radial Basis Function (RBF) approximation model prior to the identification step to reduce the computation time required to identify an optimum value of m [30]. As Fig. 8 shows, the RBF approximations are an excellent representation of the simulation results. A Downhill Simplex optimization method [30] is used to identify optimum values of the stiffness using the force response function for five independent experiments. The identified optimum values, mopt, are listed in Table 1. The average stiffness identified using the inverse method was mopt = 1.05 ± 0.18 (k = 2). Thus, within the observed repeatability level, the material is not strain rate sensitive, as m = 1.0 represents the quasi-static response. The 17 % uncertainty in mopt reflects random errors due to experimental and modeling approximation factors.
Fig. 7 Effect of friction coefficient, f, and stiffness, m, on (a) 0F and (b) 0S
Fig. 8 Typical response surface showing the effect of m on actual 0F and actual 0S along with RBF approximations used for calculating the optimal m
Table 1. Identified values of mopt
|
Experiment |
Average Strain Rate [s-1] |
 |
 |
 |
 |
|
1786 |
340 |
1.18 |
104 |
0.00623 |
6.33 |
|
1791 |
397 |
1.04 |
475 |
0.00607 |
6.55 |
|
1818 |
405 |
1.09 |
649 |
0.00573 |
6.38 |
|
1819 |
429 |
0.94 |
1346 |
0.00686 |
8.20 |
|
1820 |
422 |
1.02 |
897 |
0.00713 |
8.05 |
|
Average |
399 |
1.05 |
694 |
0.0064 |
7.10 |
|
U* |
70 |
0.18 |
930 |
0.0012 |
1.86 |
SUMMARY
An inverse method was used to determine the dynamic stiffness of a prospective biomimetic elastomer using a Kolsky bar with the intention to determine whether the strain rate sensitivity mimics real tissues. High-speed digital image correlation (DIC) was used to measure the surface deformation of the sample during the test, which suffered from significant friction effects and non-uniaxial deformation due to the extreme softness of the elastomer. Further, a special flyer technique was used to compress the specimen to larger compressive strains that would otherwise have been possible for the NIST Kolsky bar. A finite element model of the experiment was constructed and the resulting simulations were compared to the DIC shape history data and to force history data. The sensitivity of the model to specimen stiffness, friction, Poisson ratio and damping were examined, and the simulation results were compared to the experimental data to identify the dynamic stiffness using an inverse method. An appropriate numerical damping coefficient was determined that produced realistic-looking force-history signals while avoiding adding artificial stiffness to the simulated specimen which would confound attempts to identify the true strain-energy-dependent dynamic stiffness. A Poisson ratio of 0.495 was chosen to achieve reasonable computation times, but this led to a systematic under-prediction of the outer surface displacement of the sample compared to the DIC measurements. Friction in the model was very important for capturing the actual deformation of the specimen. Interestingly, the sensitivity of the simulation to friction above f = 0.5 was minimal, so a no-slip friction condition was used. With the friction and damping conditions established, the stiffness scale factor m was observed to produce distinct minima in the force and shape residual plots in five repeat experiments. The optimal stiffness scale factor for these high strain rate tests was m = 1.05 ± 0.18, indicating very little rate sensitivity in this material up to 60 % compressive strain. Thus this particular elastomer is not strain rate sensitive, unlike what has been reported for actual soft tissues.
