## Abstract

The problem is considered of a fibre that is driven dynamically, by compression at one end, into a matrix. The fibre is not initially bonded to the matrix, so that its motion is resisted solely by friction. Prior work based on simplified models has shown that the combination of inertial effects and friction acting over long domains of the fibre–matrix interface gives rise to behaviour that is far more complex than in the well-known static loading problem. The front velocity may depart significantly from the bar wave speed and regimes of slip, slip/stick and reverse slip can exist for different material choices and loading rates. Here more realistic numerical simulations and detailed observations of dynamic displacement fields in a model push-in experiment are used to seek more complete understanding of the problem. The prior results are at least partly confirmed, especially the ability of simple shear-lag theory to predict front velocities and gross features of the deformation. Some other fundamental aspects are newly revealed, including oscillations in the interface stresses during loading; and suggestions of unstable, possibly chaotic interface conditions during unloading. Consideration of the experiments and two different orders of model suggest that the tentatively characterized chaotic phenomena may arise because of the essential nonlinearity of friction, that the shear traction changes discontinuously with the sense of the motion, rather than being associated with the details of the constitutive law that is assumed for the friction. This contrasts with recent understanding of instability and ill-posedness at interfaces loaded uniformly in time, where the nature of the assumed friction law dominates the outcome.

## 1. Introduction

This paper will address certain aspects of the more general problem illustrated in figure 1. In the general problem, a fibre that is initially bonded to a matrix is loaded dynamically in either tension or compression at its end, with a load that is some function of time. The fibre and matrix are elastically dissimilar. The assembly may be of plane or axisymmetric geometry. Boundary conditions for the axisymmetric case may be those of type I or type II (i.e. stress or displacement conditions at the boundary of a cylindrical cell), as designated by Hutchinson & Jensen (1990) to simulate an array of fibres; or those corresponding to an isolated fibre in a semi-infinite body. Boundary conditions for a plane problem may be periodic to simulate an array; or again those of an isolated fibre in a semi-infinite body. A state of initial compression may exist across the fibre–matrix interface, e.g. due to residual thermal stresses. The load causes the fibre to debond progressively from the matrix, allowing relative sliding between the two, which is opposed by friction. The constitutive law of the friction may involve displacement, displacement rate and the magnitude of the interface compression. The debonding mechanism may result in a nonlinear process zone in the debond crack wake of significant length. Friction will generally act over much longer lengths, often the whole domain of sliding.

In this paper, model experiments and numerical simulations of the push-in problem (compressive end loading) will be used to study the particular limit that the interface is initially debonded (limit of zero debond toughness). One motivation is to explore the connections between the physics of long range frictional sliding and shear fracture, using model experiments in which a reduced set of mechanisms operates. Simultaneously, the question is addressed of how simply dynamic interface failure in composites can be represented, without loss of accuracy. The latter query is inspired by seminal experiments in the static case on cylindrical fibres in a matrix (e.g. Marshall 1992). In the static case, shear-lag analysis of the fibre/matrix debonding problem, which assumes Lamé-like field solutions, proves very accurate for most parameter ranges (Hutchinson & Jensen 1990) and forms an immediate link to the problem of composite fracture (Marshall *et al*. 1985; Budiansky *et al*. 1986). Where static shear-lag models are accurate, analytical results usually stand independently of whether the fibres are more or less compliant than the matrix and whether plane or axisymmetric conditions exist; preliminary shear-lag analyses of the dynamic case suggest the same compass, although the physics are considerably more complicated than in the static case (Sridhar *et al*. 2003). The problem of figure 1 may also be considered as a particular case of delamination failure in a symmetrically laid up laminated composite, for which the plane geometry is directly applicable and the lamina represented by the fibre may be more or less compliant than the other layers. The questions of how accurate shear-lag theory might be in dynamic cases, and what subset of the physics of interface failure might remain reasonably well represented in it, are therefore of practical interest. Issues specific to the push-in problem will be addressed here, including how the degrees of freedom used in representing displacement fields in modelling and the assumed nature of the friction law influence stability, front formation, front velocity, stress fields, etc.

### (a) Friction acting along bars

Shear-lag models of the dynamic case of the problem of figure 1, with friction that is spatially and temporally uniform in magnitude (but not in sign), have revealed the regimes of behaviour that might be expected for an initially debonded interface (Nikitin & Tyurekhodgaev 1990; Cox *et al*. 2001; Sridhar *et al*. 2003). First, the front of the furthest deformation in the fibre does not travel at the longitudinal wave speed for the fibre, *C*_{f}, but at some other velocity that depends on the friction strength, the elastic mismatch, the wave speeds in the fibre and the matrix and the loading rate. The front velocity is bounded by *C*_{f} and the longitudinal wave speed in the matrix, *C*_{m}, and may be lower or higher than *C*_{f}. Second, depending on the same parameters, distinct regimes of slip, slip and stick and reverse slip exist, even when the load point displacement is monotonic. Reverse slip refers to the condition that the relative motion of the fibre and matrix is opposite in sense to that of the load point. Reverse slip is possible when the wave speed in the matrix is higher than that in the fibre, so that the fibre is pulled along, at the deformation front, by the matrix. The domains of slip, stick and reverse slip are fixed (growing similarly with time) for linearly increasing loads, but appear and disappear in complicated sequences for other loading cases, even those that are simple functions of time (e.g. step functions). The solutions are much more complicated than for the same problem in static loading, where the only history dependence that affects the solution is the order of any load reversals (e.g. Marshall 1992). There is even a hint (unproven) that the birth and death of slip, reverse slip and stick-slip domains may be chaotic, with arbitrarily small changes in the load point history resulting, at later times, in finitely dissimilar domain patterns. Numerical work has supported the shear-lag results for linearly increasing loading, under the modelling conditions assumed, for most parameter values (Sridhar *et al*. 2003). However, the strong nonlinearity of the friction law, which may reverse sense at moving boundaries, makes accurate numerical location of the boundaries very challenging, even under the constraint of assumed Lamé-like solutions (N. Sridhar, K. Dunn & B. N. Cox 2002, unpublished work).

### (b) Friction between half-spaces

The problem of an interface with no debond energy but long, possibly infinite, domains of friction has been studied extensively in the geometry of two remotely and uniformly loaded half-spaces, contacting on a plane. Among a number of insightful articles, the introductory sections of Cochard & Rice (2000) provide an especially helpful summary. Some salient results are as follows. (i) For elastically dissimilar materials subject to Coulomb friction, uniform slip becomes unstable and slip pulses can form instead, i.e. slip domains on either side of which the interface is not sliding (Weertman 1980; Adams 1998). Slip pulses can propagate in either direction, in different cases, and at values of the remote shear stress that are less than the friction stress. (ii) Unstable slip is also possible for elastically homogeneous cases, but not for a simple Coulomb law. (iii) Cochard & Rice go on to distinguish instability from ill-posedness: the former exists when spatial perturbations to the state of slip grow in amplitude; the latter refers to the non-existence of solutions of any kind. Instability does not preclude the possibility of uniform slip as a formal solution, given the hypothetical presence of perfectly uniform conditions. Ill-posedness points to an inconsistency in the physics of the problem as stated. (iv) Whether the interface problem exhibits ill-posedness or instability depends on both the degree of elastic mismatch and the nature of the assumed friction law. The relevant measure of elastic mismatch is the point at which the generalized Rayleigh velocity (a mode of interface wave propagation found for a frictionless interface along which loss of contact does not occur) ceases to be defined. Ill-posedness prevails for all friction coefficients when the elastic mismatch is mild (generalized Rayleigh wave exists); and for friction coefficients above a modest critical value when the mismatch is more severe. Thus, ill-posedness is the common case. (v) Ranjith & Rice (2001) assign blame to the simplified physics of Coulomb's law, in particular its implication that friction can change discontinuously in time upon a finite change in the normal contact pressure. A modified law, in which the change occurs smoothly over a characteristic response time, had in fact already been inferred from experiments (Prakash & Clifton 1993; Prakash 1998). Such a law regularizes the ill-posed planar interface problem; and the regularized solutions continue to exhibit slip-pulse character (Cochard & Rice 2000).

The complexity seen in the Lamé-like solutions to the push-in problem of figure 1 has different physical origins to the instability and ill-posedness discussed in the last two paragraphs. In the work based on Lamé-like solutions, the friction force was assumed to be uniform and constant, apart from sign reversals, and not influenced by changes in the normal compression at the interface. Complexity arises from the effects of friction on the propagation of wave pulses along the fibre and matrix system, excited by the dynamic end loading, and is associated with the birth and death of slip, stick and reverse-slip domains, which is complicated by the strength of the nonlinear effects due to the sign-reversal property of the friction. In the problem of planar frictional interface between half-spaces, instability and ill-posedness are direct consequences of the assumed relation between friction and the normal traction, and depend strongly on whether Coulomb's or another law is assumed.

The final point from the body of work on planar interfaces that is of present relevance is the relation of the boundedness of slip-pulse domains to the nature of the assumed friction law. Zheng & Rice (1998) looked at this question with a velocity-weakening friction law, expressed as a decrease in the coefficient of friction as the shear displacement rate rises. (This law is similar to that of Prakash & Clifton, but expressed in inverse form, in terms of the displacement rate rather than the rate of change of the shear traction.) For such a law, slip motion can occur either as a crack-like event, in which the interface slips simultaneously over a semi-infinite domain extending back from a slip front; or as a slip-pulse motion, with slip confined to finite, moving domains. Which case prevails (in the problem of contacting half-spaces) depends on the rate of velocity weakening: if the slope of the decline is small, crack-like slip behaviour occurs; if it is large, slip-pulse motion occurs (Zheng & Rice 1998).

### (c) Effects of interface bonding

A separate, major body of literature has been written on the problem of interface shear crack propagation in systems in which the interface is initially bonded. The most pertinent results here refer to the question of allowed and preferred crack velocities. Energy considerations suggest that, unlike mode I cracks, which are forbidden above the Rayleigh wave speed, shear cracks can in principle propagate at any velocity, *c*, in the so-called intersonic regime, *C*_{S}<*c*<*C*_{L}, where *C*_{S} and *C*_{L} are the shear and longitudinal wave speeds (Broberg 1994; Burridge *et al*. 1979). However, numerical studies imply a restrictive condition for intersonic propagation, that the crack tip process zone be diffused over a finite interval and not concentrated at a point, as in a classical brittle crack (Andrews 1976). In the work of Burridge *et al*., the velocity √2*C*_{S} assumes a special role: for *c*<√2*C*_{S}, crack acceleration is unstable (associated with decreasing load); while for *c*>√2*C*_{S}, crack acceleration is stable (increasing load).

These early theoretical deductions have now been complemented by experiments that have achieved shear crack conditions for the first time in the laboratory by exploiting the propensity of laminated materials to confine cracks to prescribed planes (Rosakis *et al*. 1999). Specimens consisting of monolithic slabs laminated at a bonded interface as well as composite ply laminates have been studied. Under appropriate dynamic loading, approximately mode II conditions can be maintained, in contrast to tests with a uniform body (free of pre-existing weak planes), where crack deflection or bifurcation would quickly cause a transition to mode I conditions. The laminate experiments confirm mode II crack propagation in the intersonic regime and are also consistent with the predicted stability transition at √2*C*_{S}, since the crack is observed to accelerate rapidly past the shear wave speed, but decelerate and achieve approximately uniform velocity at √2*C*_{S}. Needleman has performed numerical simulations of this particular experimental configuration, where the driving force is a dynamic impact end-load rather than the remote shear loading of the early theoretical studies, and confirmed most aspects of the observed behaviour (Needleman 1999). Needleman's model incorporates mixed mode cohesive elements to represent the debonding process, with cohesive laws (traction versus displacement) that feature both an initial hardening phase and a decaying tail. These laws are somewhat different from that in Burridge *et al*., which was a purely softening law and for mode II displacements only. The details of the law do affect the outcome in some ways. For example, while Needleman's simulations show the same qualitative history of rapid acceleration to a nearly constant velocity as seen in experiments (Rosakis *et al*. 1999), the value of the final velocity generally exceeds √2*C*_{S}, rather than equalling it, with the speed attained depending on the maximum stress assumed for the cohesive law. In a separate theoretical study of shear crack propagation along bi-material (Homalite/steel) interfaces, a large jump in terminal crack velocity was found at a certain prescribed load point velocity (Needleman & Rosakis 1999). For inferior load point velocities, crack propagation was limited by the Rayleigh wave speed in the Homalite (the slower medium); above the transition, the terminal velocity approximately doubled, to be near the longitudinal wave speed of the Homalite.

While intersonic propagation is physically admissible, shear crack propagation remains forbidden in the velocity interval, *C*_{R}<*c*<*C*_{S}, where *C*_{R} is the Rayleigh velocity. The mechanism by which a crack can jump across this velocity gap involves the creation of a daughter crack ahead of the parent crack, which then links to the parent in an effective surge of crack growth (Burridge 1973; Andrews 1976; Gao *et al*. 2001). While relationship of this phenomenon to slip-pulse solutions for planar frictional interfaces must exist, the presence of a velocity gap must be associated with the presence of concentrations of tractions where interface debonding is taking place, which do not generally appear in interfaces coupled solely by friction.

### (d) Contrasts in physics

Thus, shear crack propagation analysis, pulse-slip solutions for the planar frictional interface, and Lamé-like solutions of the push-in problem with zero debond energy (figure 1) all show that variable front velocities up to the maximum of the longitudinal wave speeds in the fibre and the matrix are physically admissible. However, important distinctions arise between the different problems. In particular, there is no analogue in the Lamé-like solutions of the push-in problem with zero debond energy to the velocity gap in the interval (*C*_{R}, *C*_{S}) in fracture problems, nor any special role for the velocity √2*C*_{S}. Neither do direct analogues of the phenomena of instability and ill-posedness that have been so extensively studied for planar frictional interfaces appear in the existing Lamé-like solutions, since these have been developed for uniform and constant (not Coulomb) friction.

Solutions in the literature for dynamic shear cracks and the initially-debonded push-in problem also differ in another very important way. Even though cohesive models have been used for dynamic shear crack modelling, thus spreading the debond process into a diffuse process zone, the dimensions of this process or cohesive zone have been relatively small. For example, in the simulations of Needleman (Needleman 1999), the cohesive zone is approximately 1 mm (a result of the assumed parameter values in the cohesive law), while the crack propagates over distances of approximately 25 mm. Behind the cohesive zone, the fracture surfaces remain traction free. A small process zone representing debonding is very different from the conditions expected in the presence of friction, especially in mode II problems, where shear tractions may act over the whole cracked specimen. While a shear cohesive traction enters the fracture problem in exactly the same way that a frictional traction would enter, the extension of cohesive modelling to friction zones that may be as long as the crack has not been made in dynamic fracture studies. (Calculations for bi-material interfaces did predict the existence of large-scale contact zones trailing the crack-tip (Needleman & Rosakis 1999), but friction in these zones was not modelled.) Models have either treated long zones of friction in the absence of a debond process (the Lamé-like push-in problem or the planar frictional interface problem), or debonding as a relatively short process zone in the absence of long zones of friction.

Final introductory comments are directed upon the nature of the cohesive law. In static fracture of a material that is not rate dependent, only two or three degrees of freedom in the cohesive law have a significant influence on the fracture response (e.g. Bao & Suo 1992; Cox & Marshall 1994; Massabò & Cox 1999). If the zone has finite extent, most characteristics of crack propagation are determined by two parameters, which may be taken as the critical traction in the cohesive law and the area under the law (work of fracture). If the cohesive law incorporates an initial hardening phase (traction rising with displacement), then the cohesive mechanism may act over the whole crack, even for very long cracks, in which case a third parameter, the slope of the hardening phase, controls crack growth (Cox & Marshall 1994). If the cohesive mechanism exhibits rate dependence (excluding inertial effects), e.g. due to creep or viscosity, then one or at most two additional parameters, describing the rate of decay of the cohesive tractions, appear to suffice to predict fracture behaviour (Cox & Sridhar 2001). The possibilities when inertial effects enter may be richer. For example, in none of the work on static fracture are there cases where the rate of decline of the law with increasing displacement plays a significant role (except as it affects the area under the curve, i.e. the work of fracture). Yet in dynamic studies of friction, the rate of decline of the coefficient of friction, i.e. of the cohesive shear traction due to friction, with displacement rate has a profound effect: it determines whether slip occurs in a crack-like or slip-pulse mode (Zheng & Rice 1998). This and other results, including theoretical models (Lapusta & Rice 2003) and novel laboratory simulations of slip-pulse motion (Xia *et al*. 2004), show that constitutive laws for dynamic friction need to be represented with more degrees of freedom than laws in static facture. Inertial effects appear to manifest further details of the traction law in fracture behaviour.

All of these background works have addressed some aspects of the general problem of figure 1, but they remain separate pieces of understanding, not yet connected to one another. Thus, large-domain friction problems have not yet embraced the influence of a strong bond (energetically significant crack tip process); while fracture work has not yet extended to very large (possibly semi-infinite) domains of friction. Generality in the friction law has begun to appear only in studies devoted to friction as a phenomenon, but these studies imply that over-simplifying friction laws can have a strong affect on the possible behaviour in other systems. And last, the most general studies of fibre push-in as a special problem, with general end-loading histories, remain to be undertaken.

## 2. Experiments

A model planar specimen geometry was designed to investigate dynamic fibre push-in for a bi-material system in which no chemical bond exists at the material interface. A transparent and birefringent polymer sheet (Homalite-100), acting as a ‘fibre,’ was placed between two steel plates, analogous to a ‘matrix’ (figure 2). The dynamic fibre push-in process was mimicked by loading the fibre with a projectile travelling at speeds ranging from 10 to 40 m s^{−1}. The resistance of the polymer/metal interfaces to slip was manipulated by applying a static compressive pre-stress () acting perpendicular to the interfaces, whose value varied from 0.5 to 175 MPa.

Discussion of the resulting deformation will refer to the coordinate system, (*x*_{1}, *x*_{2}), with the origin located at the mid-plane of the Homalite piece, at the end where it is impacted (figures 1 and 2). The half-height of the Homalite will be denoted *h*. Thus, the interfaces lie at *x*_{2}=±*h*. The origin of time, *t*, is the moment of first impact.

*Materials and specimens*. The densities and relevant wave speeds of the test materials are given in table 1. All plates were 9.5 mm thick and 82.5 mm in length. The steel pieces were 50.8 mm in height. The height, 2*h*, of the Homalite piece was 16.5 mm. The Homalite dimensions were chosen to provide a fibre aspect ratio, defined as length divided by the half-height, *h*, of 10. The steel height was chosen to minimize the interactions of wave reflections from the specimen boundaries with the dynamic sliding process. This assures that the case studied is equivalent to that of an isolated fibre in a semi-infinite body. The pieces of the specimen were aligned and stacked vertically (steel/Homalite/steel) without the use of adhesive or bonding agent.

*Dynamic loading and characterization*. The fibre was loaded dynamically by impacting its end with a 25.4 mm diameter cylindrical steel projectile, 44.5 mm in length, which was accelerated using a light air gun. Projectile speeds were varied over the range 10–40 m s^{−1} by varying the gun pressure. A small steel tab was affixed to one end of the Homalite piece to eliminate the possibility of impact damage to the end of the relatively brittle polymer during dynamic loading (figure 2). The dynamic sliding of the Homalite relative to the steel was characterized using high-speed photography in conjunction with dynamic photoelasticity. A strain gauge was mounted to the surface of the steel tab to provide a trigger signal for the initiation of the recording of high-speed photographs (described below). From sequences of images, the dynamic loading history, *σ*(*t*), of the fibre, the initiation time for the onset of interface sliding and the interface slip speed histories were obtained.

Dynamic photoelasticity apparatus was configured as follows. A collimated laser beam, 100 mm in diameter, was passed through a circular polarizer, the length of the Homalite, and a second polarizer to generate the photoelastic images, which were recorded using a high-speed digital camera (Cordin 220-16, Cordin Scientific Imaging, Salt Lake City, UT 84119, USA). The camera was capable of recording 16 frames. The frames were recorded every few microseconds; the interframe times were varied in each experiment. In dynamic photoelasticity, the fringes show constant-value contours of the difference of principal stresses, *σ*_{1}−*σ*_{2}.

### (a) General observations

A total of 16 experiments were conducted. The role of impact speed (loading rate) was investigated at a constant static pre-load; and that of pre-load was investigated at a constant impact speed. Figure 3 shows a representative high-speed isochromatic fringe pattern recorded during sliding in the common member of these two test sets (, *v*_{0}=38±2 m s^{−1}). The approximately vertical fringes represent the stress generated from the impact (from the left) and propagate to the right. Shortly behind the propagating front, kinks form in the fringes in angled bands. These bands intersect the top and bottom interfaces, implying non-smooth variations in the interface displacement. Among the many fringe patterns collected from all tests, consistency in the pattern of bands was imperfect. In some images, either noise or experimental variance made it difficult to confirm or refute the presence of bands that might have been expected to be present based on their presence in other images. In the images with better-defined fringes, two bands were prominent, typified by the dashed white lines in figure 3. Based on changes in the spacing or angle of fringes before and after a band, or in the derivative of the spacing, the first of these is interpreted as the beginning of a domain of increasing interface shear stress; the second as the point of transition to interfacial sliding. In figure 3, the domain marked ‘increasing friction’ (or interfacial shear stress) is clearly defined as such by the fringe behaviour; in the domain marked ‘constant friction,’ the fringe variations are consistent with this interpretation, but the implication is not as clear due to noise. A peak in the fringes on the specimen mid-plane locates the onset of unloading.

Figure 4 presents images similar to that in figure 3 but taken from experiments conducted with higher and lower static pre-loads. The difference in the visibility of bands in these figures and figure 3 is representative. Further inspection of figures 3 and 4 shows that as the pre-load increases, the kinking of the fringes is associated with larger distortions, suggesting that the onset of sliding is associated with a sharper discontinuity. A similar trend was observed as the impact speed was increased. The dynamic sliding process is associated with higher, more localized stresses at higher loading rates and higher static pre-load.

### (b) Sliding speeds

The velocity of the dominant band in each image along the interface was estimated by comparing its position in successive frames. Figure 5 shows plots of the positions of the intersections of a symmetric pair of such bands with the top and bottom interfaces as a function of time for an experiment conducted with 12.5 MPa pre-load and 38 m s^{−1} impact speed. The bands propagate at the same velocity on the top and bottom and the variation of sliding position with time is essentially linear. This contrasts with the common case for shear crack propagation, where a period of acceleration of the crack tip is often observed, which may approach constant velocity asymptotically. Here, the velocity of the bands is constant throughout the experiment, within the experimental error of approximately ±3% (typically on the order of ±50 m s^{−1}). These observations were consistent in all of the experiments, indicating that on this scale of observation, sliding is a steady-state process. In this particular experiment, the slope of the plot in figure 5 indicates a sliding speed of 2050 m s^{−1} or approximately the dilatational wave speed of the Homalite. Model analysis (see below) suggests that these bands probably corresponded to the beginning of a domain of enhanced interfacial shear stress.

Figure 6 provides a summary of the velocities of the dominant bands measured from all of the experiments, plotted as a function of static pre-stress. While the data show substantial scatter, individual points are accurately determined (see above). In all cases, the observed speeds were in excess of the shear wave speed of the Homalite, i.e. greater than 60% of the dilatational wave speed. In general, the trend is toward decreasing sliding speed with increasing static pre-load. However, there is no clear functional dependence. Insight provided by modelling (see below) suggests that the bands identified as dominant in collating the data of figure 6 may have corresponded in some cases to the onset of interfacial sliding, but in other cases to the beginning of a domain of increasing interfacial shear stress.

Data for all of the experiments conducted at 0.5 MPa and varying impact speed lie along the ordinate (the axis of sliding speed) in figure 6. There was no systematic variation of sliding speed with impact speed, with all the measured sliding speeds being approximately the dilatational wave speed of the Homalite.

### (c) Sliding stresses and loading rate

The photoelastic fringes can be equated to the local difference in principal stresses(2.1)where *F*_{σ} is the stress-optic coefficient (22.6 kN m^{−1} for the Homalite), *N* is the fringe order, and *W* is the specimen width (in the through-thickness direction, *x*_{3}). The stress along the centre line of the Homalite can be determined from equation (2.1) alone if the following conditions are met: (i) the principal stresses along the centre line act in the *x*_{1} and *x*_{2} directions and; (ii) plane stress conditions exist; and (iii) strains *ϵ*_{22} in the *x*_{2} direction are negligible. Then equation (2.1) reduces to(2.2)where *ν* is Poisson's ratio for the Homalite, and hence the stress distribution along the centre of the fibre can be determined. The condition of negligible strain component, *ϵ*_{22}, will not be met in the presence of large compressive pre-stress. Nevertheless, equation (2.2) proves useful, for the following reasons. Numerical simulations (see below) show that, in the absence of compressive pre-stress, dynamic contributions to *ϵ*_{22} are very small and the stress, *σ*_{1}, estimated from experiments using equation (2.2) and that calculated by the simulations are very close. When the compressive pre-stress is large, its contribution to *ϵ*_{22} will be invariant in space and time. Therefore, equation (2.2) can be used to estimate the applied loading rate, d*σ*_{1}/d*x*_{1}, in steady-state conditions from experimental data, since this involves taking a derivative.

The variation of stress with position is found to be approximately linear, yielding a constant value of d*σ*_{1}/d*x*_{1} for each test. (Thus, experiments and simulations (see below) confirm the rapid attainment of steady-state conditions in the axial stress and the proper application of equation (2.2).) This steady-state value can be divided by the dilatational wave speed in the Homalite to give a loading rate per unit time. The loading rate rises in proportion to the impact speed and decreases linearly with static pre-load (figure 7). When the coupling of the fibre and the matrix is increased, the input momentum is transferred more rapidly to the matrix, resulting in slower loading of the fibre.

## 3. Numerical modelling

A computational simulation of the test was set up in a commercial finite element package (ABAQUS1), with cohesive elements to represent interface friction. Symmetrical boundary conditions (*u*_{3}=0; *τ*_{12}=0) were imposed along the centre of the Homalite piece (*x*_{2}=0) and the top surface of the steel piece. A state of plane stress was assumed in the (*x*_{1}, *x*_{2}) plane. Semi-infinite elements were used to prevent elastic waves bouncing back from the right end of the model (figure 8). The cohesive elements introduce a relationship between the shear tractions, *τ*_{int}, and the sliding displacement rate, , at the interface. The law used (figure 9) is intended to represent friction that is uniform and constant in magnitude, but an initial linear part is included to avoid numerical difficulties where the sense of the motion changes. A dynamic loading history (figure 10), corresponding to that caused by the experimental impactor, was imposed as a prescribed stress, *σ*_{11}(*t*), acting uniformly over the interval −*h*<*x*_{2}<*h* at the left end of the Homalite. The simulation described in detail in the following had a loading rate, d*σ*_{11}/d*t*=31.5 MPa μs^{−1}, a rise time of 6.0 μs, and linear unloading over a further 12.0 μs. The loading rate and maximum stress correspond to those deduced from measurements of stress distributions at the end of the Homalite for the test with an impactor velocity of 39 m s^{−1} and a compressive pre-stress, . The unloading rate was set to be half that of the loading rate in magnitude, which corresponded approximately to the measured rate. Thus, the boundary stress was not predicted, but fitted to experiments; but all other characteristics of the deformation, including front velocities and the general distribution of stresses, etc. were predicted given only this fitted condition. The interfacial frictional stress was assigned the magnitude *τ*_{f}=40 MPa. This would correspond to a friction coefficient of 0.8 for the same test, if friction had obeyed Coulomb's law in the test; but this, of course, may not have been the case. The magnitude of *τ*_{f} was held constant in the simulations. Since the loading rate depends only weakly on the compressive pre-stress (figure 7), the simulation could also be interpreted as representing cases of higher compressive pre-stress, if a lower value of the coefficient was believed to be more appropriate. (The calculation would be unchanged because the coefficient of friction does not enter explicitly into the model.)

### (a) General character of the predicted wave deformation

Figure 11 shows predicted contours of the principal stress difference *σ*_{1}−*σ*_{2} and the shear stress *τ*_{12} in the Homalite at time *t*=9.822 μs. These plots and plots of stress variations along single lines (see below) distinguish four zones in the wave propagation. Rightmost is the so-called head wave zone, where the interface friction stress is opposite in sense to that in the trailing zones. (This change of sign is not evident in figure 11*b*, because of the small magnitudes of the stresses involved.) In the head wave zone, the Homalite is loaded by the steel and the wave speed exceeds that of the Homalite, *C*_{f}, by a factor of more than 2 (the factor being difficult to pinpoint numerically) but is less than that of the steel, *C*_{m}.

Following the head wave zone is a linearly increasing shear stress zone, where the interface (friction) stress increases approximately linearly with position up to the limiting value, *τ*_{f}, given by the cohesive law. This zone derives from the finite-sloped transition from negative to positive shear stress used in the cohesive model to avoid numerical problems (figure 9). However, domains with this character were also observed in the model experiments, suggesting that the law of figure 9 may fortuitously represent, at least qualitatively, the true constitutive behaviour of the interface under rising shear stress. The next zone is the constant friction zone, which is characterized by a suddenly and greatly reduced density of the contours in *σ*_{1}−*σ*_{2}, while the shear stress is constant at the value, *τ*_{f}.

The last zone is associated with the unloading part of the loading history. The contour lines near the interface become chaotic, which is also seen in the experiments (figures 3 and 4). In the simulation, alternating contact and separation zones exist along the interface in this region. The experimental contours show similar variations, suggesting the same contact behaviour (although in the experiments no direct confirmation of interface displacements is possible).

### (b) Stress distribution and front velocities

Comparison of the measured and predicted stress contours has shown that, even though the numerical model is based on an *ad hoc* and probably incorrect friction law, most aspects of the experiments appear to be reproduced. Further insight into the physics of the deformation is available from other calculated characteristics.

Figure 12*a*–*c* shows stress profiles at three different instants during the simulations of figure 11. Figure 12*a* shows the axial stress along the centre line of the Homalite, *σ*_{11}(*x*_{1}, 0), along with a measurement from the experiment executed with impactor velocity, *v*_{0}=39 m s^{−1} and compressive pre-stress, . The agreement between the simulation and the experiment is quite close. (See also the discussion following equation (2.2).) The furthest disturbance evident in this plot corresponds to the transition in figure 11 from the head wave to the domain of linearly increasing friction. The peak axial stress along the centre line decreases as the wave propagates, which is attributed to the dissipative frictional sliding process along the interface. Invariance of the stress profile confirms the attainment of steady state conditions by time 6.2 μs, shortly after peak load is attained at 6.0 μs.

Figure 12*b* shows the interface shear stress at the same three instants. The shear stress increases approximately linearly until the constant friction stress assigned in the cohesive law is reached. It remains almost constant through the ‘constant friction’ zone of figure 11, although small oscillations are evident around the prescribed value, *τ*_{f}. The unloading zone is dominated by high frequency, severe fluctuations in *τ*_{int}.

Figure 12*c* shows the axial stress, *σ*_{11}(*x*_{1}, *h*), along the interface. A knee behind the furthest visible disturbance corresponds to the transition from linearly increasing to constant interfacial stress. In the region ahead of the knee, the loading rate of the material is remarkably small.

Figure 12*d* illustrates further the nature of the interface fields in the chaotic unloading zone. Large spikes in the fibre sliding displacement periodically reach down to the sliding displacement of the matrix, indicating locations of fibre/matrix stick. Reverse slip is also possible for brief intervals at these locations, but is difficult to resolve in a numerical simulation, due to computational noise.

The plots of figure 12 provide reasonably accurate information for evaluating the velocities of the various fronts. Figure 13 shows histories of the locations of fronts as functions of time, obtained by interpolating measurements of the locations of the corresponding features in figure 12. The leading edge of the head wave is not shown, since it could not be accurately discerned.

The velocity, *V*_{2}, of the transition to the linearly increasing friction zone is evaluated from both the earliest significant axial stress along the centre line and the earliest significant shear stress along the interface. The former gives *V*_{2}=1.04*C*_{f}, while the latter yields a slightly lower value *V*_{2}=0.98*C*_{f}. To within reasonable uncertainty, neither is significantly different from the wave speed in the Homalite.

The front of the constant friction zone can only be determined accurately from the interface stresses, either the knee in the axial stress or the onset of the shear plateau (figure 12*b*). The inferred wave velocity is *V*_{3}=0.87±0.02*C*_{f}, the error reflecting the difference between the two sources and the axial stress implying the higher value.

The front associated with the onset of unloading is best identified by the peak in the axial stress along the centre line of the Homalite (figure 12*c*). The velocity takes the value *V*_{4}=0.98±0.02*C*_{f}, once again close to the wave speed in the Homalite.

The velocity estimates can be approximately confirmed by analysing the bands or kinks in the predicted stress contours, which propagate from locations on the interface at which a transition (non-smooth behaviour) is present in the displacement fields. Two bands are traced by dashed lines in figure 11, emanating approximately from the fronts associated with the onset of linearly increasing friction and constant friction (sliding). The angles subtended by the lines to the interface, *α*_{i}, *i*=2 and 3, are given by sin *α*_{i}=*C*_{S}/*V*_{i}, where *C*_{S} is the shear wave speed in the Homalite. While difficulty arises in locating the bands on some contours, the fitted values sin *α*_{2}≈0.45 and sin *α*_{3}≈0.6 are probably correct to within 20%. These values yield *V*_{2}=(2.2±0.4)*C*_{S} and *V*_{3}=(1.6±0.3)*C*_{S}, in reasonable agreement with the more accurate values deduced above from figure 12. The boundary of the domain of chaotic fields is also marked by a dashed line in figure 11, which emanates from the unloading front. The angle subtended by this line, *α*_{4}≈arcsin 0.35, yields *V*_{chaos}=(3.0±0.3)*C*_{S}. This high implied velocity is consistent with the observation in the simulations that the onset of chaos is delayed somewhat after the beginning of unloading, but the chaotic zone tends subsequently to catch up with the unloading front, and therefore propagates at a velocity exceeding that of the unloading front, which is approximately *C*_{f}.

The different velocities predicted for the onset of sliding and the beginning of the domain of increasing interface shear help understand the distribution of the data in figure 6. Some uncertainty arose in identifying bands in fringe images, so that the nature of the interfacial phenomenon corresponding to a particular band was not always clear. Given the velocities predicted by the model, those data in figure 6 falling close to *C*_{f} would be associated with the beginning of the domain of rising interfacial shear stress, while other data would be associated with the sliding front.

Since the front of the constant friction zone propagates at slightly different velocities at the interface and along the centre of the Homalite, the configuration of contours must change with time. Other contour plots at different times (not shown) show that the nearly parallel, vertical contours in the ‘linearly increasing friction’ and ‘constant friction’ zones of figure 11 do move apart more rapidly at the centre of the Homalite than at the interfaces, which causes them to bow out towards the head wave with increasing time. The same fringe divergence was seen (via dynamic photoelasticity) in the model experiments. In both the simulations and the experiments, the divergence is limited in magnitude and the stress contours remain at angles less than approximately 10° from vertical during passage of the deformation along the length of the test piece.

Other changes with time are: (i) the zone of chaotic contours grows into the interior of the Homalite with time. In fact, other contour plots from the same simulations as reported in figure 11 reveal that chaotic behaviour does not appear at all until the elapse of approximately 8 μs. (ii) The constant friction stress zone shrinks with time, since *V*_{4}>*V*_{3}, indicating that the unloading wave propagates faster than waves ahead of it.

## 4. Results from shear-lag models based on Lamé fields

For the static problem of figure 1 and related thermal and mechanical problems, many useful results have been obtained using shear-lag analysis and the assumption that the fibre deformation fields are separable in *x* and the transverse or radial variables; and that fields are simple, known functions of the latter (Marshall *et al*. 1985; McCartney 1987; Cox 1990; Cox *et al*. 1990; Hutchinson & Jensen 1990). Analysis shows that these approximations give quite accurate predictions of the end load versus the end displacement as long as the slip distance is long compared to the fibre diameter (or fibre width in a plane problem; Hutchinson & Jensen 1990). The fibre fields are well represented, and also the matrix fields, provided the volume fraction of fibres is not too small. Since these results in the static case are simple and have enlightened many aspects of the mechanics and engineering principles of designing fibre-reinforced materials, the question of their utility in the dynamic case is of interest.

For dynamic end loading that increases linearly in time and friction that is uniform and constant, shear-lag analysis predicts that the deformation will propagate in two or three domains, depending on the value taken by the following three dimensionless parameters (Sridhar *et al*. 2003):(4.1)where *f* is the fibre volume fraction. The specimen of figure 2 and the simulations have lateral boundary conditions that are not periodic and therefore are not in strict correspondence with the shear-lag model of Sridhar *et al*., but approximate correspondence can be sought, using the dimensions of the specimen, by setting *f*=8.25/50.8=0.162. The domains predicted by shear-lag modelling consist of different interfacial conditions: slip, stick and reverse slip. For the case of the simulations reported in §3, *k*=0.317, *C*=0.378 and *ϕ*=0.0043. Shear-lag modelling predicts, for these values, that the deformation during the linearly rising part of the loading will comprise a stick zone and a slip zone (see appendix A). Furthermore, the interface friction stress in the stick zone is negative and has a very small magnitude, 0.036 MPa (see equation (A 6)). The sign is negative because the matrix (steel), having the higher wave speed, is pulling the fibre along in this zone. The magnitude is small because the volume fraction of the fibre is small, so that the stresses in the matrix remain small, the spatial derivative remains small, and therefore the friction stress remains small. For the case studied, the shear-lag model predicts (appendix A) that the head wave front will propagate at *V*_{1}=2.6*C*_{f}.

A further prediction of the shear-lag model is that, if the fibre volume fraction rises, then the interface shear stress in the head wave zone will exceed the friction stress in magnitude and reverse slip will occur (consider figure 14 as the parameter *ϕ* rises). The head wave zone would no longer be a stick zone.

Figure 11 broadly confirms the predictions concerning the head wave zone. The zone is a stick zone (no slip), the interface shear stress is negative and very small in magnitude, and the front velocity, while not well determined, is consistent with the value 2.6*C*_{f}.

The zone of linearly increasing friction in figure 11 has no analogue in the shear-lag results of Sridhar *et al*., because in that particular modelling, the friction stress was assumed to change sign with reversing slip direction as a step function. The zone of linearly increasing friction in figure 11 arises from the linear friction stress regime at small displacements in the law of figure 9. The constant friction zone of figure 11 is equivalent to the slip zone in the shear-lag model. The front of the constant friction zone is predicted in the shear-lag model to propagate at 0.73*C*_{f}. This is less than the velocity, *V*_{3}=0.87*C*_{f}, computed in the simulations of figure 11. An interesting question is whether the two models would agree in the predicted front velocity for higher fibre volume fractions, equivalent lateral boundary conditions, and as the linearly varying domain in the constitutive law of figure 9 shrinks in extent.

If it is assumed that the friction stress value used in the shear-lag model would vary with compressive pre-stress according to Coulomb's law, then a prediction of the variation of front velocities with compressive pre-stress can be made. Two such predictions for the velocity of the sliding front, *V*_{3}, have been superimposed on figure 6. The rates of change of the friction stress correspond to assuming a coefficient of friction of 0.1 or 0.4. The trends of the curves are not inconsistent with the weakly defined declining trend in the data, especially when those data lying at *C*_{f}, which are believed to correspond to the front of increasing interfacial shear stress, are excluded.

## 5. Possible chaotic behaviour

The noisy behaviour appearing in the unloading zone in the simulations does not appear to be a numerical artefact, but rather a consequence of the modelling assumptions and the physics of the problem. If it were due to Gibbs' phenomenon, i.e. numerical ringing at a near-discontinuity in the solution, it would not be confined to the vicinity of the interface (note, for example, the smoothness of the axial stress at the centre line of the Homalite, figure 12*a*); nor to the unloading zone, since non-smoothness of equal strength exists in the fields at other fronts (zone boundaries).

The noisy behaviour is not associated with Coulomb's law, since this was not the constitutive assumption of the model. In fact, the friction law assumed (figure 9), which contains no relation between the friction stress and the normal pressure, is simpler than that used in most studies of either stability or ill-posedness in problems of slip between half-spaces. For the same reason, the predicted phenomenon cannot be directly related to Poisson's effect, even though Poisson's ratio was non-zero in the simulations.

The fact that the noisy behaviour is confined to the unloading zone, while the deformation propagates smoothly at a range of speeds in other zones, has no obvious analogue in any prior studies of instability. The closest to a hint of similar behaviour is in incomplete and unpublished results for the shear-lag modelling described by Sridhar *et al*. (2003). Attempts to trace solutions for linear unloading that follows linear loading lead to apparently complex patterns of the birth and death of slip, stick and reverse-slip zones, as the unloading wave overtakes deformation of earlier origin. The possibility is suggested (but not yet thoroughly researched) that arbitrarily small changes in the loading history cause finite changes in the subsequent pattern of slip, stick and reverse-slip zones. In the shear-lag problem, the only possible source of instability or complexity is the strong nonlinearity associated with the change in sign of the friction stress when the direction of interfacial slip changes. While the analysis is incomplete, the suggestion of chaotic behaviour in the shear-lag model provokes the speculation that the noisy behaviour seen during unloading in the experiments and the numerical simulations is also an example of chaos.

Shear-lag modelling also highlights the role, in the possible generation of chaotic behaviour, of the coupling of wave motions between the fibre and the matrix. For a fibre in a rigid matrix, where the fibre deformation is modified by interfacial friction tractions that depend only on its own motion, stable and well-behaved solutions can be mapped out systematically for quite complex loading and unloading histories, e.g. by the method of characteristics (Nikitin & Tyurekhodgaev 1990). This method is not useful when the fibre motion is coupled to that of a matrix.

Numerical solution of the shear-lag problem, e.g. numerical solving of difference equations subject to the constraint that displacement fields are Lamé-like, is not enlightening, because the rapid births and deaths of different zones are quickly obscured in numerical noise at the fronts.

The similarity of the predicted domains of chaotic behaviour and those measured in the model experiments is very eye-catching (figures 4 and 11), but the conditions assumed for the interface are unlikely to be the same in the two cases. In particular, the friction law of figure 9 is an idealization that is likely to be an oversimplification of the true friction behaviour in the experiments (although the experimental law is not known in detail). This suggests that the putative existence of chaotic behaviour is not especially sensitive to details of the friction law, but could be a consequence of unloading with any friction law that possesses the strong nonlinearity of sign reversal.

## 6. Conclusions

Comparing the characteristics of model planar experiments, numerical simulations and shear-lag modelling leads to the following inferences about the push-in problem of interest (figure 1), for a bi-linear sequence of loading and unloading.

The broad characteristics of the deformation in the fibre, exclusive of shock waves, can be predicted at least qualitatively by simple shear-lag models, in which the displacements are reduced to Lamé fields. The shear-lag solutions have the advantage of proposing clearly defined zones of different slip behaviour, which while present in both the experiments and the numerical simulations, are not so easily identified in them without prior knowledge. For engineering studies of dynamic deformation in composites reinforced by fibres, rods, or stitches, shear-lag results will convey at least the trends of dynamic response at the macroscopic level. For example, the dynamic fibre end displacement is likely to be described quite well by shear-lag theory, since those fine details of the push-in (or pullout) phenomenon that shear-lag theory misses are unlikely to strongly influence the end displacement, which is derived from integrated strains.

For the given geometrical conditions and loading rate, shear-lag modelling correctly predicts the system of stick followed by slip found in the numerical simulations and model experiments, even though the details of the assumed or actual friction laws are likely to be different in all three. One infers that the pattern of zones is mainly determined by the interplay of stress waves of different speeds in the fibre and the matrix and is insensitive to the friction constitutive law.

The speeds of fronts (or the boundaries between zones) are fairly well predicted by shear-lag modelling, but more consistent between the numerical simulations and the model experiments. Shear-lag modelling could be expected to be a better approximation as the volume fraction of the fibres rises.

As far as shear-lag theory is accurate, results found here for plane geometry should be equally applicable to axisymmetric conditions, which are a useful approximation to composites of cylindrical fibres.

The literature is replete with analyses of frictional sliding between half-spaces under uniform far-field velocity conditions. For that problem, the presence of ill-posedness, instability, pulse and crack-like slip systems, etc. depends strongly on the nature of the friction law (see §1). The present study considers non-uniform (time-dependent) loading, which yields results that are not accounted for by the prior literature. Complex behaviour, which may be chaotic, arises during unloading, but not at other stages of the deformation, in both the experiments and the numerical solutions. This behaviour appears to be a consequence of the strong nonlinearity associated with the sign reversal in the frictional tractions with slip direction and is not sensitive to the details of the friction constitutive law.

The results that some characteristics of engineering relevance, especially the fibre end displacement and the pattern and velocity of sliding fronts, are insensitive to details of the friction law (being so similar in the experiments, numerical modelling and shear-lag analysis) also contrasts with the literature on shear crack propagation. One would expect that as the frictional traction rises in magnitude, or a large debond energy is included in the problem, the fibre push-in problem would begin to manifest the complex dependence seen for shear cracks. At the same time, the conditions required for shear-lag analysis to be accurate would be lost.

Changing the compressive contact pressure in the experiments has quantitative but not qualitative effect on the deformation patterns.

## Acknowledgments

Q.D.Y. and B.N.C. are grateful to support from the U.S. Army Research Office through contract number DAAD19-99-C-0042, administered by Dr David Stepp.

## Footnotes

↵ABAQUS, Inc., Pawtucket, RI 02860, USA.

- Received April 22, 2005.
- Accepted October 31, 2005.

- © 2006 The Royal Society