Phase transition modeling of polytetrafluoroethylene during Taylor impact

Multi-phase and multi-component materials including polymers and polymer matrix composites have attracted the attention of researchers in the last several decades following the wide use of these materials in aerospace, defense, and automotive applications, which requires prediction of the response of such materials under extreme conditions. Polymers are well known to exhibit multi-phase behavior that is easy to observe in experiments, as well as strong strain-rate dependence. Moreover, understanding the metastable phase transition behavior is important due to similar mechanisms observed in geological material such as quartz and its derivatives. The factor that complicates analysis and needs to be taken into account is that the phases in condensed multi-phase materials usually have non-negligible strength.

Polytetrafluoroethylene (PTFE), also known under the trademark Teflon, is a model polymer material, possesses a simple atomic structure containing one amorphous and four known crystalline pressure and temperature dependent phases. At low pressures, the material transits through three low-pressure hexagonally packed crystalline phases, the atomic structure of which are distinguished by changes of ordered or disordered molecular rotations controlled by temperature. At high pressures, the low-pressure crystalline phases transform into a crystalline phase (phase III) characterized by a planar atomic structure and orthorhombic or monoclinic lattice packing. Numerous publications (e.g., Refs. 1 and 2) have reported on the high-pressure transition with observed critical pressure P_cr between 0.5 and 0.7 GPa. A few indications of this high-pressure phase transition in the plane wave impact tests were observed as Hugoniot cusps^3–5 corresponding to the transition. Moreover, more detailed micro-structural study⁶ has revealed that crystallinity increase follows the high-pressure phase transition in shock loaded material.

When experimentally studying PTFE under impact conditions, an anomalous rapid ductile-to-brittle transition has been observed^7,8 in Taylor cylinder impact tests on increasing the impact velocity, U₀, above a critical value, U_c. This effect has been interpreted as the manifestation of the high-pressure phase transition. At lower velocities, the PTFE rod exhibits the classic mushrooming compression-induced bulging at the front part of the Taylor specimen near the anvil as the compression wave propagates along the rod, as observed^7,8 in many polymers in Taylor cylinder impact tests. For most polymers, as the impact velocity increases experimental observations^9–14 find the initiation and arrest of stable ductile cracks. These stable cracks can occur at the centerline of the specimen, petalling as the tensile hoop stresses at the circumference of the front surface are large enough to initiate fracture, or as spiraling shear cracks from combined tensile hoop stresses and compressive axial loading. For most polymers, these cracks grow longer with more cracks initiating as the impact velocity increases until the cracks start linking up and the rod end fails catastrophically. This classic transition from purely ductile deformation to diffuse brittle fracture occurs over a range of multiple tens of m/s. The PTFE Taylor impact specimens on the other hand exhibit non-classical response where the transition from apparent ductile-to-brittle behavior consistently occurs over a range of impact velocity of only 1 m/s. This impact velocity U₀ = U_c of the observed ductile-to-brittle transition corresponds to the phase III transition pressure, P_cr, in PTFE as estimated⁷ with hydrocode simulation. Phenomenological effects of the phase transition in PTFE on the velocity profiles in planar impact, which are not as easy to measure in situ during the Taylor cylinder impact test ductile-to-brittle transition have been observed experimentally in Ref. 8, and constitutive modeling of a three-phase approximation of PTFE with amorphous and two crystalline phases has been conducted in Ref. 15.

The present paper approximates the PTFE polymer as a two-phase material and focuses on the most significant change of PTFE mechanical properties due to the phase III transition and increase of crystallinity. The polymer under normal conditions is a mixture of amorphous and crystalline phases. The equilibrium pressure-temperature phase diagram of PTFE is well documented (see Ref. 16), whereas the metastable phase behavior at high pressure and under rapidly changing loading and thermal conditions in shock waves (including a hysteresis in crystallinity) has yet to be determined and is being studied at present. A low-pressure crystalline phase is phase II under normal conditions migrating to phase IV and further to phase I with increasing temperature, which operates by untwisting helical molecules of the crystalline phases into less consolidated forms. For the present two-phase representation, the mechanical change in impedance between low- and high-pressure phases linked with the phase transition is of primary interest, because the main objective of this paper is the mechanical response to the impact loads in the multiaxial mode of loading found in Taylor tests. The mechanical variations that result in changes in impedance stem from changes in elastic modulus, densities, and yield limits that are associated with the phase transition. Density variations are not great between the low-pressure crystalline phases and only the rise in crystallinity has major effect upon the density. With the pressure increase, the crystalline phases transform first into a planar molecular modification of phase III, whilst the amorphous phase transforms into a crystalline one. Again, the nature of the transition from the low-pressure to the high-pressure phase does not need to be directly specified, i.e., it could be the change in mechanical characteristics due to either an increase in crystallinity or the transition into phase III since either option may result in similar mechanical properties. The only relevant information needed is the thermo-mechanical properties of the phases and a transition kinetic linking the phases. Keeping in mind that mechanically, the density of the low-pressure crystalline phases (the phases II, IV, and I) and the crystallinity are not changing significantly, we may assume that the mechanical properties of the homogenized low-pressure phase (the first phase in our approximation) are close to those of the bulk material (the polymer under normal conditions), which is represented as a mixture of the amorphous and one or several of the three low-pressure crystalline phases. Alternatively, the high-pressure phase (the second phase in our approximation) is assumed to be close to phase III since the crystallinity increases. The kinetic of the transition from the first to the second phase is mostly governed by pressure with a slight dependence of the transition zone on temperature following the phase diagram. The mechanical mismatch in impedance driving the phase transition is manifested in different mechanical properties of this phase to the material under normal conditions.

Numerical analysis employing the two-phase model implemented in the CTH hydrocode is conducted in the present paper for the tests⁷ complemented by X-ray digital tomography of PTFE samples from similar Taylor tests at U₀ = 106 m/s and U₀ = 135 m/s. The analysis links the observed brittle failure mechanism for PTFE in the Taylor tests at above the critical impact velocity with the phase transition driven by the mechanical impedance difference. At the impact velocity U₀ below the critical velocity U_c, the material may be fractured by the viscous damage mechanism typical for polymers. At the same time, studies⁶ have hypothesized that the abrupt start of the brittle fracture mechanism at U₀ = U_c is associated with a sudden volume change when density jumps from its initial value for an amorphous microstructure up to that typical of a crystalline one as the phase transition proceeds. At the same time, another study¹⁷ observes formation of the crystalline fibrils surrounded by flaws at crack tips with increasing deformation rate supporting the hypothesis of an association of the nucleation of the damage spots with the crystallinity rise in the molecular cluster domains. Following this low-pressure experimental analysis, it can be also assumed that the crystalline domains of higher packed phase III are promoting the damage at high pressure and strain rates. To describe this mechanism, a simplified kinetic of failure associated with the phase change is used for modeling the damage.

The two-phase model has been developed¹⁸ in a form suitable for hydrocode implementation. For convenient model implementation, the structure of hydrocodes, such as CTH,¹⁹ requires a certain desirable model framework. For example, the P-lambda model by Grady²⁰ implemented in the hydrocode has such a structure associated with the evolution of internal state variables along the material particle trajectory. The present two-phase model was originally formulated in Ref. 21 and employed for moderately porous materials. Further development²² of this model neglected the strength of the phases but considered the inter-phase internal energy exchange, which enabled us to use the model for highly porous metallic powders and treat adequately anomalous Hugoniot behavior. Strength was taken into account¹⁸ in a recent version of the model. A more complex three-phase version has been developed¹⁵ and used for constitutive modeling of PTFE plate impact tests that allowed us to adequately describe experimental velocity histories.⁸ The version¹⁸ of the two-phase model with included strength has been implemented^23,24 in the CTH hydrocode for porous materials. The present paper is a further implementation of this model in CTH for the case of two condensed phases linked by a transition.

The present approximation of PTFE as a two-phase material used in the CTH model implementation is employed in the present paper in order to analyze the phase transition effects on the damage and failure of PTFE rods during Taylor cylinder impact tests.

The PTFE polymer is considered as a mixture of amorphous and crystalline domains under normal conditions (phase 1 of the two-phase material representation) and a high-pressure phase approximating the planar crystalline phase III of the polymer (phase 2 of the two-phase material representation). The underlying two-phase model¹⁸ with strength contains the mass, momentum, and energy conservation laws. The conservation laws are complemented with the equation of state outlined below and the constitutive equations. The latter are for a number of internal state variables, namely, elastic deviatoric strains to describe the plastic flow, inter-phase strain disbalance, entropy disbalance, and the mass and volume concentrations.

The model equations are closed with the equations of state (EOSs) for each of the phases presented in a form that describes internal energy as a function of density, deviatoric strains, and entropy. The mass weighted internal energy, e, for the material is represented by a mixture of the phases in a total EOS formulation for the two-phase materials thus

where c is the mass concentration of the first phase, and superscripts refer a variable to the corresponding phase. Focusing on the choice of EOS for a condensed two-phase material, we specify each phase individually by choosing the energy potentials e^(k)(ρ^(k), e_ij^(k), s^(k)) (k = 1, 2) in the compact form^15,18,25 of the Mie-Gruneisen-type EOS. Here, ρ^(k), e_ij^(k), and s^(k) are density, mass weighted strain components,¹⁸ and specific entropy in the corresponding phases. The EOS constants are specified in Ref. 15. The exponent constants for the modulus are obtained from available Hugoniots. The thermal capacities c_v, Gruneisen parameters γ, initial densities ρ₀, and initial values of the bulk and shear modulus K and G will be specified below. In some cases, the moduli are determined from the phase elastic Young's modulus E₁ and E₂ and Poisson ratios, ν₁ and ν₂.

The elastic moduli in PTFE are highly temperature sensitive as observed in Ref. 26 and analyzed in Ref. 15. However, in the present Taylor tests,⁸ temperature only increases by 10°–20°; therefore, in the present analysis, moduli were selected to be constant. The corresponding Young's modulus and Poisson's ratio for phase 1 are taken as those for the PTFE in as received condition and E₁ = 0.88 GPa and ν₁ = 0.47. The initial densities of the phases are assumed to be ρ₀₁ = 2.17 g/cm³, ρ₀₂ = 2.335 g/cm³, and density of the material in the initial state, ρ₀ = 2.175 g/cm³. Similarly,¹⁵ the Young's modulus for the second phase is taken to be 30% higher than that for the phase 1, E₂ = 1.146 GPa with the same value of Poisson's ratio ν₂ = 0.47. The corresponding bulk and shear modules are determined from G = E/[2(1 + ν)], K = E/[3(1 − 2ν)]. The remaining EOS constants are also taken from Ref. 15, for phase 1 represented as weighted average of “amorphous” and low-pressure crystalline (“cs”) phases and for phase 2 as a high-pressure crystalline (“ch”) phase using representations.¹⁵ 

The constitutive equation for the phase transition kinetic is taken in the following thermodynamically consistent form:²² 

Here, d/dt is the particle trajectory derivative, Λ = μ₁ − μ₂ is an analogue of the chemical potential providing a stationary point for c when the specific Gibbs energies, μ, of the first and second phases are equal. The Gibbs energy specified for the first phase, for instance, is μ₁ = e⁽¹⁾ + p⁽¹⁾/ρ⁽¹⁾ − T⁽¹⁾s⁽¹⁾ − e_ij⁽¹⁾s_ij⁽¹⁾, where for this phase, p⁽¹⁾ is pressure, T⁽¹⁾ is temperature, s_ij⁽¹⁾ are components of the stress deviator tensor, and e_ij⁽¹⁾ are components of a tensor representing an extensive measure of elastic deviatoric deformations.¹⁸ The exchange rate for the mass concentration c is a simplified version of the function¹⁸ representing Avrami kinetics with a transition pressure, P_cr, slightly decreasing with temperature in accordance with the PTFE phase diagram¹⁶ 

where $p1′$ = p⁽¹⁾ − p′ and $p1″$ = p⁽¹⁾ − p″. Here, p′(T⁽¹⁾) is a linear fit correlating with the boundary of the low-pressure phase transition. This phase transition is actually neglected in the present work for the present two-phase representation. The associated Arrhenius term is kept in Eq. (3) only for consistency with the kinetic.¹⁵ With the present choice of kinetic, this term practically does not affect variation of the mass concentration responsible for the phase transition. This insensitivity is in agreement with the present representation of phase 1 as an average of the crystalline phases I, IV, or II mixed with the amorphous phase. As in Ref. 15, the kinetics are fixed to have two different rates for the phase transformation in the forward and reverse directions. In the forward transition, a polymorphic mechanism is assumed, whereas in the reverse transition, a martensitic mechanism is accepted (discussion on the reasons for this can be found in Ref. 15). A line p″(T⁽¹⁾) correlates with the high-pressure (phase III) phase transition,¹⁵ which is relevant to the present phase choice in contrast to the line p′. The corresponding polynomial coefficient in Eq. (3) is taken to be zero at p⁽¹⁾ < p″, and s1 is the second shear stress invariant. The exponents in Eq. (3) are n_ε = 0.01, m₀ = 3, and q₀ = 1. The activation energy of the Arrhenius term in Eq. (3), U(p) = U⁰ − f₀ × (p/ρ₀₁), has the following coefficients: U⁰ = 8.32 J/g and f₀ = 0.018. The remaining material constants in Eq. (3) are selected as follows: B₀ = 6.55 × 10³ s/m², C₀ = 44 s/m², c_min = 0.001, c_max = 0.999, p₀ = 1 atm, and R = 8.314472 J/K/g. The characteristic constants for the p′ and p″ lines are listed in Table I.

The characteristic pressure constants in Table I correspond to the p″-line that represents the high-pressure phase transition in the kinetic (2) and (3), and this is slightly lower than the values from the wide observation range of 0.45–0.79 GPa reported in the experiments.^1,4 Briefly, this apparent discrepancy can be explained by the strength of the material, which results in the difference between the corresponding component of stress recorded in experiments and pressure used for the transition in the kinetic equations. This will be discussed in detail in a subsequent section.

The compaction kinetics regulates the volume concentration of the first phase, θ, by inter-phase forces equilibrating the concentration according to the following equation:²² 

Here, Π₀ is the potential¹⁸ tending to equilibrate phase pressures. The rate for the parameter θ is evaluated from the material compressibility similarly to the previous choice¹⁸ when simulating the α-ε phase transition in iron

Here, p is pressure in the two-phase material, A₀ = 0.01, n₀ = 1.05, and p_c = 3 GPa.

The compaction kinetics (4) describes the change of the phase volume concentrations due to mechanical loading of material with invariable phase content. The phase transition introduces an extra factor in the volume change, which must be taken into account. The phase transition contribution, ψ_ph, which is an addition to the mechanical compaction rate, ψ, in Eq. (4), and which is applicable to the mass and volume concentrations c and θ below 1 (including the concentrations close to zero) has been derived in Ref. 18. However, in the present case the concentrations can be close to both 0 and 1. Therefore, we need to generalize the correction ψ_ph for the arbitrary case.

Let us consider mechanical deformation of a representative volume with mass m followed by phase transition. We consider the process during a time increment Δt. The total mass is unchanged, however, the total volume may change due to both deformation and phase transition in each of the phases, resulting in the total rate change in Eq. (4) as −ψ + ψ_ph. Thus, the volume V = V₁ + V₂, containing the mass m = m₁ + m₂ (the subscripts refer to phases 1 and 2) is incremented by ΔV = ΔV₁ + ΔV₂ with denotation for the incremented values V′ = V + ΔV, so that V′_k = V_k + ΔV_k (k = 1,2) with incremented masses m′₁ = m₁ + Δm and m′₂ = m₂ − Δm.

Now, we can separate the change of volume due to deformation that evolves according to Eq. (4) from the change due to phase transition. In the analysis that follows further, we consider the volume changed only due to the phase transition. As a result, density in the phases does not change because the mass increment in first phase Δm (Δm can be either positive or negative) corresponds exactly to the volume increment ΔV₁ with similar logic for the second phase. Using the definitions of mass and volume concentrations during this process, we have

where

Thus, the original, V_k, and incremented volumes, V′_k, along with the total one, V′, can be calculated from Eq. (6) via masses and phase densities. Using those, the volume concentration increment is

where Δc = Δm/m. After some manipulations with the density definitions (6), the volume concentration increment reduces to

For this process related solely to phase transition, dividing the volume concentration increment by Δt with Δt converging to 0 gives the corresponding contribution ψ_ph to the volume concentration rate

Combination of the compaction rate from Eq. (4) with the phase transition contribution transforms Eq. (4) to

with ψ₀ from Eq. (5).

The heat exchange rate for relevant exchange parameter (inter-phase entropy disbalance) is obtained from evaluation of heat equilibration between phase clusters.^18,22 The heat transfer coefficient, h, for this rate is taken from Ref. 18 

Here, d₀ is a characteristic dimension, which is of the order of phase domain size, and k₁ and k₂ are the thermal conductivities of the phases. A_S is the dimensionless Nusselt number associated with the phase domain factors such as hydrodynamic convection within a phase. These factors may promote or resist inter-phase heat exchange. As discussed,¹⁵ A_S is of the order of 1 in the present case. The size of crystalline fibrils in the damage zone associated with the brittle fracture can be in a quite wide range between tens of nanometers and hundreds of micrometers.^17,27 Choice of this parameter will be discussed in detail in the modeling section. As a preliminary baseline for the calculations, we select a lower bound d₀ = 10⁻⁸ m.

The yield flow is governed by the stress relaxation functions. The functions are determined²⁸ from two yield stresses (for example, the dynamic, Y_d, and static, Y_s, yield stresses) versus two strain rates for each of the phases. The values of Y_s and Y_d at normal (T = 20 °C) and elevated (T = 50 °C) temperatures^26,29 for phases 1 and 2 adopted in this case are listed in Table II.

The damage mechanism is assumed to be associated with the volume change by introduction of flaws and voids when the crystalline fibrils are formed at the tip of crack as analyzed in Ref. 17. From this assumption, the damage can be characterized by a variable ξ as introduced in the following equation:

Here, ξ is the damage parameter that manages the introduction of voids into a material according to the approach used in the CTH hydrocode.¹⁹ The function D in Eq. (9) is chosen in the simplest form D(c) = 0.01 × H(c − c_m)s⁻¹, where c_m = 0.5, and H is the Heaviside step function. It should be noted that thermodynamics of mechanical and internal variables is not subject to the change of ξ. The only effect of the parameter increasing from zero to one is an increment of the failure volume in the CTH algorithm resulting in the voids growing within the material.

The model implementation in the CTH hydrocode¹⁹ follows similar steps²³ to those described for two-phase porous material model. The input parameters are initialized in the Lagrangian module of the code. During the time cycle calculations, the basic and the internal state variables responsible for the inter-phase exchange processes (mass and volume concentrations, entropy disbalance, strains, etc.) are accessed from the Lagrangian and Eulerian modules. The calculation starts in the Lagrangian module with the deviatoric stress calculations for strain and strain imbalance parameters and continues in the same block for mass and volume concentrations according to Eqs. (2), (7), and (9). The heat exchange calculations using Eq. (8) start at the launch of the Eulerian module. The equation of state is a mass weighted combination of the Mie-Gruneisen-type EOSs from Eq. (1) calculated in the main Eulerian block of the CTH hydrocode.

The CTH implementation employing the kinetic functions defined above will be used in the calculations presented below.

The preliminary stage of the study is a detailed analysis of wave propagation from anvil to PTFE in the one-dimensional planar impact geometry. We consider the reverse impact by a thick steel flyer plate with velocity U₀ against a PTFE target.

If the target were a conventional single-phase material, the shock transition diagram would simply be a single Hugoniot. However, for materials subject to internal kinetic processes, the shock transition involves multiple Hugoniots with quasi-static changes of internal kinetic variables on each of them. For example, the shock transition between Hugoniots has been analyzed previously in porous materials in quasi-equilibrium²² and in Teflon represented as a two-phase mixture of amorphous and crystalline phases.¹⁵ In the present case, the representative phases have similar mechanical properties, therefore, the difference in phase equilibration is not as noticeable on the Hugoniot plane as for the representation involving the amorphous phase.¹⁵ However, for the two-phase material we can employ the mixture rule combining Hugoniots of individual phases. Schematically, the individual phase Hugoniots, H₁ and H₂, are drawn in Fig. 1 with a schematic connection path T between the Hugoniots via the phase transition process, which is detailed in the vicinity of the transition area in the inset of Fig. 1.

The phase transition connecting the Hugoniots H₁ and H₂ is rate sensitive and the curve T describing the transition is governed by the kinetic equations (2) and (3). We consider several cases corresponding to various conditions of loading with the final states indicated by points 1–5 in Fig. 1. The initial state is marked by a bold point at V = V₀₁. The impact velocity values close to U_c are selected below for the one-dimensional numerical analysis and they correspond to those used in the Taylor cylinder impact tests⁸ recorded with high-speed photography and to the present tests analyzed with X-ray tomography.

We start with the impact velocity U₀ = 106 m/s below the critical value, U_c. The final state behind the shock front for this case corresponds to point 1 in the Hugoniot diagram in Fig. 1. The calculated results are shown in Fig. 2 and demonstrate that a phase transition does not occur in this case and the wave has conventional structure with mass concentration for the first case (homogenized PTFE at low pressure) close to 1 throughout the whole wave propagation process. Density in the main wave does not reach the reference density of the second phase, and pressure is below the critical value, P_cr, initiating the phase transition. Vertical dashed lines indicate boundaries of the moving steel flyer plate and the PTFE target sample. Coordinates of the lines in Fig. 2 are displaced by the corresponding interface velocities integrated over the time difference for two different instances. These lines diverge from one another more noticeably at greater impact velocity, as will be observed in the illustrations below.

With increase of the impact velocity up to U₀ = 125 m/s, the critical pressure is just achieved, which corresponds to point 2 in the diagram of Fig. 1. However, the time to complete the transition is not enough due to a low transformation rate given by Eq. (2). As seen in Fig. 3, at t = 1 μs the transition is not yet complete, although it does so before t = 3 μs. Observing the development of the stress, it is seen that above the critical pressure within the phase 1 initiates the transition followed by a relaxation that results from equilibrating the phase pressures to a lower equal value along the Hugoniot H₂. This gives the peaks followed by relaxation in density and stress profiles at t = 3 and 5 μs as seen in Fig. 3. This relaxation gives rise to tensile stresses after reflection of the shock wave in the flyer plate from the anvil's free surface, which may be seen at t = 5 μs.

When increasing the impact velocity further up to U₀ = 140 m/s, we reach the final state at point 3 in the schematic of Fig. 1. The pressure level is still within the transition zone as seen in the calculation results in Fig. 4. However, the transition occurs fast enough that the pressure relaxation effects are barely noticeable showing a significantly higher rate of the phase transition than in the previous case. The shock is strong enough that even after the inter-phase pressure relaxation, the pressure in the two-phase material follows the Hugoniot H₁, the transition curve T, and the Hugoniot H₂ monotonically.

If the velocity is increased up to U₀ = 200 m/s, the transition zone becomes very small and the shock wave propagating in the two-phase sample appears as an almost classical two-step wave shown in Fig. 5 with the first wave following the shock transition H₁ and the final state corresponding to point 4 on H₂ (Fig. 1). The calculation in Fig. 5 also shows pressure profiles in order to assess the role of pressure in the rate (3) driving the phase transition. It is seen that the pressure developed in the PTFE sample is notably higher than that in the steel flyer plate. In fact, continuity through the contact interface takes place only for velocity and stress in direction normal to the contact interface, whereas a pressure jump takes place for materials with differing yield stresses in contact. Thus, the yield limit is significantly higher in the steel plate than in PTFE, resulting in the observed pressure difference.

As a consequence, the observed pressure is comparable with P_cr and results in a corresponding stress that is higher by an increment proportional to the yield limit. Thus, the value of critical stress corresponding to the pressure P_cr is notably higher than P_cr. This explains the choice of P_cr that is apparently lower than the transition stress (usually referred to as pressure) observed in the experiments. Choosing stress instead of pressure for the phase transition kinetics is difficult because of the need of a directionally invariant parameter characterizing load. The role of directionality on the transition is not well understood yet but the difficulties can be seen here. At the same time, the kinetic equations are intended for use in the multi-dimensional calculations required for the Taylor cylinder impact test. In order to make the choice of a characteristic force initiating the transitions that is directionally invariant, pressure was selected as a universal variable suitable for both one- and multi-dimensional calculations with a reduced value, P_cr, corresponding to the relevant transition stress.

Finally, increasing the impact velocity to U₀ = 500 m/s, the shock wave degenerates into a single step seen in Fig. 6. In this case, similarly to the overtaking of elastic precursor by strong plastic wave in the elasto-plastic-transition, the Rayleigh line connecting the initial point 0 with the final state 5 is aligned with the Hugoniot's slopes for H₁ and H₂ (see Fig. 1). It should be noted that slow evolution of the states (stress and density) behind the shock front has the same causes as those seen in the analysis of shock compression of porous materials.^18,22 These result from inter-phase heat exchange that eventuates in a relatively slow temperature relaxation.

Comparing the mass (phase 1) concentration profiles with the gradient characterizing the phase boundary in the above calculations, the non-linear locus of the boundary velocity can be seen. With respect to the particle velocity, this speed is zero up to U₀ = 106 m/s, the quickly increasing ratio of the velocity of the final (phase transition) wave to the velocity of the precursor (the Hugoniot H₁) wave from 0.2 at U₀ = 125 m/s, through 0.4 at U₀ = 140 m/s to 0.7 at U₀ = 200 m/s and reaching 1 at U₀ = 500 m/s. Thus, in the vicinity of the critical impact velocity, a stagnation zone can be localized with a small velocity of the phase transition boundary when comparing with the main wave providing mechanical deformation of the remaining material.

The results of the present one-dimensional analysis will be used in the subsequent sections for understanding the influence of the phase transition effects on the fracture behavior of the PTFE rods in the Taylor cylinder impact tests.

The PTFE polymer has been experimentally studied in the Taylor cylinder impact test.^7,8 The key finding of the tests⁷ conducted at various impact velocities, U₀, was the revealing of a critical velocity of impact U_c that is between 133 m/s and 134 m/s for a slightly different geometry of the tests.⁸ For the impact velocity U₀ below U_c, the rod was deformed and damaged following a ductile mechanism, and above U_c fracture of the rod is brittle. The present test set-up agrees with Ref. 8 and employs PTFE rods with the length, L, to diameter, D, ratio of 5 with diameter D = 9.19 mm. In addition to the high-speed observations and data on the recovered rods,⁸ the present analysis employs experimental data on recovered rods for the tests in the same set-up at U₀ = 106 m/s and U₀ = 135 m/s.

The damage within the samples resulting from the Taylor cylinder impact tests was characterised in 3-D using X-ray microtomography. A Nikon XTH225 tomography system was used, equipped with a 225 kV (225 W) high-energy microfocus X-ray source (3 μm spot size) and a 1900 × 1500 pixel² Varian 2520 14-bit 250 mm × 200 mm flat panel detector (∼130 μm pixel pitch). An X-ray energy of 50 kV and current 180 μA gave an X-ray transmission through the sample of ∼45%. An aluminium filter of thickness 0.5 mm was placed just in front of the X-ray source and acted to remove low energy photons from the X-ray beam, thus reducing beam hardening artifacts. Due to the cone beam X-ray source geometry of the tomography system, the distance between the source and sample determined the magnification and thus the voxel size in the reconstruction. The sample was positioned relative to the source and detector such that its full length (or height, as it was positioned during scanning) was being captured in the projections. This gave a voxel size in the reconstructions of ∼30 × 30 × 30 μm³. Through a rotation of 360°, 1001 unique projections were collected with an exposure time of 1 s per projection. The projections were reconstructed into 3-D tomographic volumes using a filtered back projection algorithm. Data analysis was performed using the image processing and visualization package Avizo, used to create virtual 3-D isosurface renderings. Views of the impact interfaces of the recovered rods reconstructed by three-dimensional X-ray computed tomography are shown in Fig. 7. A very limited damage area (confined fracture in terms of the analysis⁷) is seen for the case U₀ = 106 m/s (Fig. 7(a)), and an extensive brittle damage accompanied by petalling and focused in the axial region at a depth from the impact surface can be seen at U₀ = 135 m/s (Fig. 7(b)).

Hydrocode calculations⁷ associated with the experiments have evaluated pressure developed in PTFE rods at impact with the velocity above U_c. Because the calculated pressure exceeded P_cr, the authors⁷ suggested that the change in damage mechanism is caused by possible phase transition. In the present study, the Taylor cylinder impact tests have been modeled in a two-dimensional simulation using the two-phase model implemented in CTH. The set-up assigns a coordinate X to the radial direction with the symmetry axis X = 0 and Y in the direction of impact, thus, allocating the space L > Y > 0, X < D/2 to the rod and Y ≤ 0 to the anvil. In the initial stage of the impact, the pressure is focused at the axis X = 0 near the rod-anvil interface, Y = 0. Pressure in PTFE, seen as the dark gray zones of pressure fringes in Figs. 8(a) and 8(b), exceeds P_cr at a depth from the rod-anvil interface for the case of impact velocity U₀ = 135 m/s (Fig. 8(b)) that is above the critical value. For the case where impact was below the critical impact velocity, U₀ = 106 m/s (Fig. 8(a)), pressure sufficient for the transition is reached only in a small zone on the rod axis for a short period of time followed by quick lateral release. The present pressure calculation agrees with that obtained in Ref. 7, which indicates an adequate choice of material parameters for the simulated polymer. In contrast to analysis,⁷ the phase transition can be traced directly in the present calculation. As follows from the one-dimensional analysis, the phase boundary propagates slowly at the impact velocity U₀ close to U_c and accelerates when U₀ rises. As also discussed in Sec. III, due to the higher yield flow stress of the anvil material, pressure developed in the PTFE rod is apparently higher than that in the anvil, while the normal stress is continuous across the interface. This pressure jump through the rod-anvil interface can be seen in Figs. 8(a) and 8(b).

In contrast to the one-dimensional analysis, the lateral release waves in the Taylor cylinder impact tests promote stagnation of the phase boundary due to the failing pressure. Therefore, a few microseconds after the initial stage shown in Fig. 8, the pressure in the rod-anvil interface region is mainly below P_cr and the phase interface is practically stationary, similar to the results of the one-dimensional calculations in the range of nearly critical impact velocities.

The numerical set-up of Figs. 8(a) and 8(b) corresponds to the assumption of sliding rod-anvil interface, which has been attempted in the experimental set-up⁸ where the interface was lubricated with molybdenum disulfide in order to minimize the interface friction. The calculations in Figs. 8(a) and 8(b) presume that the interface between the Teflon and the steel has no friction with the lubricant simulated as a thin layer of anvil material without strength. However, increasing pressure with impact velocity may change the lubricant's viscosity, thus, affecting the shear transfer between rod and anvil through the grease. Therefore, in order to address possible friction effect, another set of calculations has been conducted with a “sticky” (non-sliding) interface condition (normal transfer of velocity and stresses) for an analysis of the possible influence of friction on the rod shape. Comparison between the calculation for the set-ups with sliding (Fig. 8(c)) and sticky (Fig. 8(d)) interface conditions at below the critical impact velocity demonstrates the effect of friction manifested by different dissipation of the shear work on the interface indicated by the temperature increase. It is seen that in the case of sliding interface (Fig. 8(c)), the shearing is focused in the rod near the interface in the vicinity of the axis. On the contrary, in the case of the sticky condition corresponding to the infinite friction (Fig. 8(d)), the shearing mainly occurs at the rod edges due to the Poisson buckling of the rod with the rod interface movement restricted by the anvil. For cases above the critical impact velocity, similar comparison of calculated dissipation characterized by temperature elevation is shown in Figs. 8(e) and 8(f). This comparison demonstrates that the sliding condition (Fig. 8(e)) promotes dissipation of the internal energy over the contact interface due to the shearing in contrast to the sticky interface case (Fig. 8(f)) directing the dissipation upwards from the interface. These differences will be attended below when analyzing the phase boundary and failure development in the rod's material.

The interface friction is not the only uncertainty of the experiment. Another important one is the phase domain size characterized by parameter d₀ mentioned in the Sec. II. Variability of this parameter involves possible variations of the material microstructure. As mentioned previously, the lower bound of the characteristic size for the baseline case was selected as d₀ = 10 nm. As a result, the inter-phase heat transfer occurs rather quickly. However, the same experimental observations demonstrate that the crystalline domain dimensions may vary significantly throughout the deformation area. As analyzed,¹⁵ the heat transfer exchange may notably affect the phenomenological response of multi-phase polymer. Avoiding detailed analysis of the influence of the loading conditions on the domain size, we conduct a comparative upper bound numerical analysis with the phase domain size larger than that in the previous case by three orders of magnitude, d₀ = 10 μm, which is also within the range observed in the micrographs.^6,17,27

Summarizing, we model the Taylor cylinder impact test using four representative simulations with parameters listed in Table III. These simulations are based on two contact interface conditions mentioned above and two different parameters of the phase domain of the modeled material affecting the rate of internal inter-phase heat exchange. These factors, distinguishing the simulation of Table III, are associated with material and test uncertainties inescapable in the experiments. Calculated rod shapes when employing the above-listed set-ups outline major possible variations and will be discussed and compared with the available experiments below.

The CTH calculations for the four set-ups listed in Table III show the phase boundary and accumulated failure at a time moment soon after the bounce event, t = 300 μs, in Figs. 9 and 10 at below and above the critical velocities of impact, respectively. It is seen that at below the critical impact velocity, U₀ = 106 m/s (Fig. 9), the phase transition occurs only at the symmetry axis in the vicinity of the rod-anvil interface for the set-ups #1, and #3 (corresponding plots in Fig. 9). When it occurs, the phase transition is initiated at the very beginning of the process as shown in Fig. 8 and is preserved through the whole duration of the process as seen in Figs. 9 and 10. For set-up #2, the material does not undergo a phase transition (Fig. 9).

The illustrations demonstrate that development of the phase transition zone at U₀ = 106 m/s (Fig. 9) is practically absent and occurs at U₀ = 135 m/s with stagnation at some depth due to proximity of U₀ to U_c as seen in Fig. 10. The distinction of the phase transition development between the cases of sliding and sticky interface conditions can be understood by referring to the analysis of energy dissipation shown in Figs. 8(c)–8(f) for these two cases. Higher energy dissipation at the symmetry axis accompanied by a sufficiently high pressure results in the phase transition even at U₀ = 106 m/s. It should be noted that similar plane impact condition (Fig. 2) does not manifest any phase transition due to smaller temperature increase in the absence of the shearing. At the same time, the effect of the sticky interface results in the temperature increase only at the rod edges where pressure is below the critical one, P_cr. The low pressure, in spite of the high dissipation, prevents the material from a phase transition at the rod edge to occur. At the symmetry axis, initially the physical state is closer to that observed in the plane impact (Fig. 2) and the phase transition is not observed. At above the critical impact velocity, U₀ = 135 m/s, the phase transition occurs regardless of the contact condition as seen in Fig. 10. The dissipation for both contacting conditions is mainly concentrated at the rod axis. Therefore, the observed above-discussed differences in the dissipation caused by the interface conditions shown in Figs. 8(d) and 8(e) correlate with the size of the phase transition zone along the rod axis. In addition, the sticky condition promotes the plane impact state in the vicinity of the axis up to duration of approximately 4.6 mm/c ≈ 3.3 μs, where c = (K/ρ)^0.5 ≈ 1.4 mm/μs is the sound velocity in Teflon at normal conditions. As seen from Figs. 3 and 4, for the present material representation, this calculated time is sufficient for propagation of the phase boundary along the impact direction for 3–5 mm. According to Fig. 10, this phase zone penetration occurs for the sticky (infinite friction) case, whereas for the sliding case the shearing at the symmetry axis followed by the smearing off the dissipation zone quickly destroys the plane impact conditions, which reduces the phase zone propagation along the impact axis.

The failure parameter “F,” the boundary for which at F = 1, is shown in Figs. 9 and 10, representing the accumulating of voids according to the CTH algorithm and is a result of a certain number reaching the fracture criterion ξ ≥ 1 in Eq. (4). In fact, this failure parameter governs material erosion representing fracture in this Eulerian hydrocode. Such defined failure is consistent with a mass concentration change at the phase transition. Nevertheless, the selected mechanism of damage is associated only with the phase transition and the failure zone is significantly wider and occupies the whole width of the rod for the case of U₀ = 135 m/s, whereas it is concentrated only at the rod axis for the case of U₀ = 106 m/s. In general, for all the set-ups the failure zone associated with the phase transition is, if it exists, localized at the symmetry axis at below the critical impact velocity and it expands over the whole rod's contact interface at above the critical impact velocity. These calculated results agree with the X-ray tomography analysis in Fig. 7.

Due to the essentially two-dimensional nature of the numerical domain, a few hundred microseconds after impact, the lateral expansion of the rod at U₀ = 135 m/s ceases as seen in Fig. 10 because the petalling observed in the experiments (Fig. 7(b)) is an essentially three-dimensional phenomenon. In turn, the sticky condition restricts the lateral deformation even more for both below and above the critical impact velocity due to contact friction limiting the lateral expansion. At the same time, if the contacting interface is under the sliding condition, the energy dissipated as lateral deformation, specifically for above the critical impact velocity case in Fig. 10, is released as fracture observed as the voids that appear due to the criterion (9). The petalling occurring in the experiments creates surface energy during the fracture process, which slightly delays bouncing compared with experiments at approximately t = 250 μs and observed in the calculations at approximately t = 200 μs.

Experimental high-speed photography observations⁸ at U₀ = 125 m/s and U₀ = 140 m/s enable us to compare the CTH calculations using the present two-phase model with the experiment. The results are shown in Fig. 11 (U₀ = 125 m/s) and Fig. 12 (U₀ = 140 m/s) within set-up #3. The observations⁸ consist of a set of frames separated by the 24 μs time increments. The impact event starts at t = 0, and the before-impact time interval, t₋, is determined from the distance between the rod and anvil S₋ and velocity of the rod, U₀, as t₋ = S₋/U₀. Then, the time after the impact for the second frame is calculated as t₁ = 24 μs – t₋, and so on. Thus, calculated times for Fig. 11 start from t = −11 μs with the interval of 24 μs (for example, the last frame in Fig. 11 corresponds to t = 253 μs). Similarly, the starting time in Fig. 12 is t = −23 μs with the last frame corresponding to t = 241 μs.

The comparison shows agreement with the effects at moving edges, which, as mentioned before, are associated with the friction or petalling due to failure. The friction effects in set-up #4, when compared with the present calculation, are similar to those in Figs. 9 and 10. Note, friction effects manifest themselves more clearly at impact velocities below critical levels that may in fact be associated with the grease viscosity reducing at a higher contact pressures for impact velocities above the critical value. As mentioned previously, there is a restriction in edge expansion due to hoop stresses in the two-dimensional set-up in the rod-anvil contact zone at U₀ = 140 m/s caused by the petalling in experiments and this results in a smaller calculated expansion of the rod. With this exception, the results of modeling are close to experimental observations.

Comparison of the recovered rod lengths with those calculated is not straightforward. The difficulty of the hydrocode simulation of the Taylor cylinder impact test is that the rod after bouncing from the anvil continues to change its length for some time due to interaction of waves that slowly attenuate between the ends of the rod. For the case of impact with U₀ = 106 m/s, the calculated coordinates of the ends of the rod are drawn as hexagonal points at every 20 μs during the first 1 ms after impact as shown in Figs. 13 and 14 for set-up #1. Results for the remaining simulations from Table III are similar in behavior for the rod. For illustration, the corresponding results for the set-up accounting for friction (set-up #2) are shown in Figs. 13 and 14 as insets for the first half millisecond. The only significant difference between the numerical results for these set-ups is velocity of the rod as it is released. This is seen when comparing set-ups #1 and #2 in Figs. 13 and 14. The solid curve traces the length of the rod (coordinate of the top points minus coordinate of the bottom points). It is seen that the change in length follows wave equilibration along the rod corresponding to half of the period T′ with velocity c = (K/ρ)^0.5 ≈ 1.4 along the rod length L ≈ 46 mm so that T′ ≈ 2 L/c ≈ 65 μs, This agrees with the half period time interval in Figs. 13 and 14.

The start of bounce phase is approximately the same for all the set-ups at approximately t = 200 μs. In general, since the in-coming release is a result of interaction between the rod and anvil rarefaction waves, the bounce event and interfacial releases may depend on the anvil design and environment of the experiments, and, as mentioned earlier, petalling may delay the start of the bounce. The rod lengths are very similar for all the set-ups and shorten as a function of the impact velocity. The observed recovered shapes are shown in Figs. 13 and 14 to the right from the X-ray digital tomography borrowed from Fig. 7 (Fig. 13) and from photographs of recovered samples⁸ (Fig. 14). It is seen that the rod lengths are in a good agreement with the average line drawn through the middle of the oscillating curves.

In order to compare calculated rod profiles with the experimental data, we take calculations at time equivalent to those taken from the experiment. The comparison of the sliding condition calculation is drawn on the left from the axis X = 0 and the sticky condition calculation on the right from the axis for the both pairs of impact velocities; experimental (dashed) profiles are shown in Fig. 15 for the whole set of simulations from Table III.

The sliding condition, when comparing cases (#1 and #3) with (#2 and #4) in Fig. 15, does not introduce a significant shape and dimension effect for the rod, although it somewhat affects the phase transition and damage zones due to differences in release effects. Thus, the comparison of observed shapes of the recovered Taylor cylinder impact samples with the calculated results (Fig. 15) shows a reasonable agreement given the uncertainties in the mechanical properties of PTFE and, specifically, hypothesized mechanical properties of the high-pressure phase, phase III. Additionally, the discrepancy seen at higher velocity is consistent with viscoelastic relaxation observed in recovered polymer Taylor cylinder impact test specimens relative to in situ high-speed photographic observations.^9,30 The discrepancy in rod expansion is consistent with mass that would be conserved in the absence of petalling, which was neglected within the two-dimensional calculation in the present analysis.

Using a two-phase representation of semi-crystalline PTFE, the present work has analyzed crystalline phase transformation in PTFE rods during the Taylor cylinder impact test with the help of a two-phase model developed earlier and implemented in CTH. The calculated phase transition boundary correlates well with the experimentally observed failure zone in the PTFE rods confirming that the high-pressure transition may induce failure due to the volume change between the low- and high-pressure (phase III) phases.

The analysis confirms that pressure in the PTFE rods subject to the Taylor cylinder impact test with an impact velocity in the vicinity of the critical velocity U_c is close to the phase transition threshold identified in the one-dimensional plate impact numerical analysis. The numerical analysis of planar waves in modeled PTFE shows that the shock propagation is followed by a phase boundary, stagnating to an almost stationary phase interface at the impact velocity close to U_c. Anomalously narrow range of impact velocities resulting in the ductile-to-brittle fracture transition agrees well with the exponential kinetic nature of the phase transition development.

Considering experimental uncertainties, the calculations have demonstrated that the contacting conditions, specifically, those promoting energy dissipation at the rod-anvil interface, may slightly affect the shape of the recovered rods and the damage and phase boundary zones as may be expected. The calculations demonstrated slight inter-phase heat exchange effects managed by a characteristic phase dimension on the rod profiles. Comparison of calculated and experimental rod profiles has established the preference of the characteristic phase domain using the present model within the micrometer range.

Elastic reverberations continue for a long time after the rebound of the rod from the anvil, which also results in uncertainty in calculations of states within the recovered rod. Its reduced diameter at the impact face in the brittle fracture regime is a consequence of the continuum treatment of the PTFE material in the calculation, which prevents the rod end from flattening during fracture.

More data are needed in order to quantitatively better specify the phase transformation characteristics due to the shock waves. This would also help in confirming that a rise in crystallinity explains the observed features or whether this is only due to a transformation of crystalline phases. To do this, the mechanical properties of the amorphous phase and specific crystalline phases in shock need to be better diagnosed.

It is the case that the stress state is fully three-dimensional and significant changes in the shear component in the front will trigger the phase transformation in a more complex manner than present theories allow. The analysis shows further that the three-dimensional effects associated with petalling could affect the shape of the recovered rods at above the critical impact velocity and, thus a three-dimensional modeling with a specific fracture model accounting for the failure due to hoop stresses might be needed.

This work was conducted over a number of laboratories. A.D.R. is grateful to the CTH development team, Sandia National Laboratories, for their help in release of the hydrocode to DSTO. Thanks are due to the facility staff at Los Alamos National Laboratory, University of Manchester, Cranfield University, and AWE for experimental assistance and their contributions to the work. E.N.B. was supported by the Science Campaigns at Los Alamos National Laboratory operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. DE-AC52-06NA25396.