## Abstract

Relaxation testing is a fundamental tool for mechanical characterization of viscoelastic materials. Inertial effects are usually neglected when analysing these tests. However, relaxation tests involve sudden stretching of specimens, which causes propagation of waves whose effects may be significant. We study wave motion in a nonlinear elastic model specimen and derive expressions for the conditions under which loading may be considered to be quasi-static. Additionally, we derive expressions for wave properties such as wave speed and the time needed to reach a steady-state wave pattern. These expressions can be used to deduce nonlinear elastic material properties from dynamic experiments.

## 1. Introduction

### (a) Motivation

Relaxation testing is a fundamental tool for the mechanical characterization of materials that exhibit viscoelastic behaviour. In these nominally one-dimensional (1D) experiments, a specimen is stretched rapidly to a prescribed level, and measurements of the time-decaying force needed to maintain this stretch are interpreted to derive mechanical properties. Normally, this analysis involves an assumption of uniform axial strain in the specimen. However, inertial effects from the rapid displacement of the end of a specimen required by these tests give rise to wave motion. When wave amplitudes are significant, the strain distribution is non-uniform, and measured forces oscillate.

The study is motivated by a need to derive the stretch-dependent nonlinear elastic response of pure collagen from relaxation tests (Cohen *et al*. 1976; Pryse *et al*. 2003). Central to the interpretation of these data is an understanding of how waves propagate, how they affect measured forces and how their motion relates to mechanical properties.

This work studies these issues through a solution for wave motion in a nonlinear elastic, model relaxation specimen, and describes how the results can be applied to characterization of some viscoelastic materials. Although linear elastic and linear viscoelastic material models applied to the same specimens result in linear equations that can be solved with well-established methods (Achenbach 1973; Eringen & Suhubi 1975; Suf & Farris 1994; Ginsberg 2001; Barclay 2004), the nonlinear material model used in this paper results in a hyperbolic wave equation (Jeffrey & Taniuti 1964; Nayfeh & Mook 1979; Moodie *et al*. 1991; Bukiet *et al*. 1996; Thoo & Hunter 2003). Nonlinear wave equations can be solved using numerical techniques (Armen & Garnet 1976; Shamardan 1990; Bukiet *et al*. 1996) and several analytical schemes including the perturbation method (Giacaglia 1972; Skorokhod *et al*. 2002), the method of multiple scales (Kakutani *et al*. 1968; Achenbach 1973; Whitham 1974; Nayfeh & Mook 1979) and the method of characteristics (Fox 1955; Jeffrey & Taniuti 1964). The primary tools we employ to study the hyperbolic wave equation are the method of characteristics and finite-difference analysis.

### (b) Problem definition

The 1D model problem studied in this paper is a prismatic specimen of length *l*_{o} that is fixed at one end (*x*=0) and rapidly displaced an amount *u*_{o} at the other (figure 1). Of interest is the displacement *u*(*x*, *t*) as a function of time, *t*, and initial position, *x*.

#### (i) Governing equations

The momentum equation and small strain–displacement relations are(1.1)and(1.2)where *ρ* is density, *s=s*(*x*, *t*) is stress and *ϵ=ϵ*(*x*, *t*) is engineering strain. We adopt a ‘polynomial nonlinear’ elastic constitutive law for *ϵ*>0 (Van Den Abeele 1996):(1.3)where *E*_{o} and *k* are constant to be determined for a specimen. When *k*=0, linear elasticity is recovered.

Combining (1.1) and (1.3) results in a second-order quasi-linear governing equation:(1.4)We assume that , meaning that the equation is hyperbolic, and verify this assumption using the final solution.

We make equation (1.4) dimensionless by defining the following variables:where , and is the linear elastic wave speed. This results in the following governing equation:(1.5)where is dimensionless strain, and *μ*=0 for the linear case.

The dimensionless form of the constitutive law is(1.6)where represents dimensionless stress, and *μ* shows the effect of nonlinearity.

#### (ii) Boundary and initial conditions

We approximate the step displacement *u*_{o} by a linear ramp over a time *t*_{o}, so that the boundary conditions, when dimensionless time *T*>0, arewhere(1.7)in which . This results in a continuous displacement and a discontinuous velocity at the moving end of the specimen, as opposed to the discontinuous displacement and singular velocity that a true step would cause.

The initial condition is that the specimen is initially unstretched and at rest:

### (c) Overview

We use the method of characteristics (Jeffrey & Taniuti 1964) to analyse the nonlinear wave equation (equation (1.5)). To aid in interpretation of the results, we first review the technique by applying it to solve the linear wave problem in §2. When the solution is graphed in the space–time domain, it shows several distinct regimes of behaviour. The nonlinear wave problem is analysed in §3, using the method of characteristics. The nonlinear wave behaviour is more complicated, and several new divisions of the space–time domain graph emerge. Closed-form expressions are derived for the total wave energy and approximate expressions are derived for the nonlinear wave frequency and the time required to approach steady-state behaviour. We validate the approximate expressions against finite-difference predictions (§§4 and 5), then use these expressions to describe how best to design relaxation tests, and how to interpret nonlinear elastic material properties from a specimen's dynamic response (§6).

## 2. Background: method of characteristics and the linear solution

### (a) Elastic solution

This section reviews the method of characteristics by applying it to solve equation (1.5) for the case of *μ*=0. Although the authors have not found this specific solution scheme applied to this specific problem elsewhere in the literature, the linear solution is presented to aid in the interpretation of the nonlinear solution, and is not a focus of this paper.

The canonical form of equation (1.5) involves two first-order equations of two function variables (Nayfeh & Mook 1979):(2.1)and(2.2)where is a dimensionless velocity.

Equation (2.1) involves differentiation in the direction defined by and equation (2.2) involves differentiation in another direction defined by . Therefore, (2.1) and (2.2) can be written asalong ‘C−’ characteristic curves, defined by , andalong ‘C+’ characteristic curves, defined by . Integrating and noting that *c* is constant in the linear case (*c*=1), we have(2.3)and(2.4)where J− and J+ are ‘Riemann invariants’ which are constant along the C− and C+ characteristic curves, respectively.

### (b) Graphical representation

Figure 2 shows the elastic solution in a two-dimensional space–time domain. The C− and C+ curves are the dashed lines with slopes of −1 and +1. The ramp stretch at *X*=1 (line AB) sends a wave through the specimen. Energy is confined to the ‘disturbed’ regions bounded by the darkened characteristic bands, emanating from points A and B.

The solution procedure involves solving (2.3) and (2.4) for J− and J+ along each characteristic by applying the boundary conditions and initial conditions:

Along AB,

*U*=*T*/*T*_{o},*V*=1/*T*_{o}.Along OA,

*U*=0,*V=e*=0.Along BT,

*U*=1,*V=*0.Along OT,

*U*=0,*V*=0.

Three different regions are distinguishable in the graph. For all loading ramps, these regions reappear periodically in the space–time domain. *V* and *e* can be found in each of these regions, and *U*(*T*, *X* ) can be found from *V* and *e* (see appendix A). The regions are as follows.

#### (i) Undisturbed regions

Triangles like OAC, EBF and DHI represent the undisturbed state of *U*. For all characteristics in these regions, *U* is constant and J+=J−=0. For example, along the *X*-axis, *U*=0 for all *X* and *V*=0 for *X*<1, meaning that J+=J−=0 along all characteristics emanating from the *X*-axis for *X*<1. These regions contain no energy, since *V*=*e*=0.

#### (ii) Constant state regions

Quadrilaterals such as ABEC, DEFH and HGJI represent the disturbed state of *U* in which the wave energy is non-zero, and *V* and *e* are constant. For regions like ABEC, J+=0, since the C+ lines all pass through an undisturbed region. Similarly, for regions such as DEFB, in which the wave moves back towards *X*=1, J−=0 on the C− lines. In ABEC, the boundary condition on AB yields J−=2/*T*_{o}, meaning that *V*=*e*=1/*T*_{o}; in DEFB, *V*=−*e*=−1/*T*_{o}.

In the case of an ideal step, the disturbed regions collapse to lines. On these borders, *U* changes discontinuously from 0 to 1 or *vice versa*, meaning *V* and *e* are both infinite.

#### (iii) Boundary regions

We call constant state regions CED, FGH and IJK ‘boundary regions’. In these regions, waves exchange kinetic energy for potential energy. Here, *V*=0, and *e*=2/*T*_{o}.

## 3. Analysis

We now apply the method of characteristics to solve equation (1.5) for *μ*≠1 and compare the solution to the elastic case. An interesting feature of these nonlinear waves is that wave energy spreads throughout the entire domain after a time *T*_{f}. We find approximations for *T*_{f}, the wave speed and the wave frequency. Additionally, we find a closed-form expression for the total wave energy and present integral forms for the kinetic and strain energies.

### (a) Nonlinear elastic solution

We begin by transforming equation (1.5) into two first-order equations in terms of *V* and *e*:(3.1)and(3.2)As in the linear case, (3.1) and (3.2) can be cast into canonical form, along C− curves, on which d*X*/d*T*=−*c*(*e*),(3.3)and along C+ curves, on which d*X*/d*T*=+*c*(*e*),(3.4)Unlike the linear case, *c*(*e*) is not necessarily constant along the characteristics. As a consequence, the characteristics C− and C+ may be curved.

Integrating (3.3) and (3.4) yields(3.5)and(3.6)where(3.7)

### (b) Graphical representation

Figure 3 depicts the first cycle of wave motion in the space–time domain. Wave energy is again confined to a few ‘disturbed’ regions.

The domain of the disturbed regions is markedly more complicated than in the linear case. The characteristic curves, which are shown below to be lines, are not parallel (note that only a few representative characteristic curves are plotted in figure 3). This has two consequences. First, waves broaden through a ‘simple wave’ region, as discussed in §3*b*(iii). Second, some disturbed regions of the space–time domain are bounded by curves other than C+ and C− curves, as discussed below in §3*b*(i).

#### (i) Boundary conditions and nonparallel characteristics

As in the linear case, *V*, *e* (hence *ψ*), J+ and J− are all zero on OA in figure 3, so J+=0 and J−=0 on any C+ and C− curve intersecting OA. As before, any C+ curve intersecting AB also intersects OA or OC, meaning that J+=0 on these characteristics. Since *V=*1/*T*_{o} on AB, equation (3.6) yields *ψ=*1/*T*_{o} and J−*=*2/*T*_{o} on AB, and on any C− curve emanating from AB.

Solving for *e* in equation (3.9) and substituting for *ψ*, we see that the strain on AB, *e=e*_{initial}, is constant along AB:(3.8)The constant *e*_{initial} means that *c*(*e*) is constant. As a result, the C− and C+ curves are straight lines with a slope of 1/*c**, where(3.9)Equation (3.9) reduces to *c***=*1 for the linear case (*μ*=0).

However, since strain is not continuous at points A and B, the characteristic lines emanating from outside AB may have a different initial strain, and thus a different *c*(*e*). Consequently, characteristics emanating from outside AB intersect those emanating from line AB. Since strain is unique at any point, characteristics may intersect only if the strain is discontinuous, and the characteristic curves end at the intersection point. We call the locus of such intersection points ‘borderlines’ to emphasize that they are not characteristics on the space–time plot.

One example of a borderline is shown in figure 3. The lower group of C− curves (long, thin, dashed lines), emanating from OA, have a slope of −1/*c*(0)=−1; the upper group of C− curves (short, thin, dashed lines), emanating from AB, have a shallower slope of −1/*c** (>−1). *e* and *V* are discontinuous over the borderline AC. The borderline AC separates what will later be shown to be an undisturbed region (OAC) and a disturbed (constant state) region (ABE_{1}C). AC also represents the wavefront in the first cycle. As will be discussed later, AC must lie between the line AC″, drawn from A with a slope of −1/*c**, and the line AC′, drawn from A with a slope of −1.

#### (ii) Wave speed

The slope of the borderline AC represents the speed of the wavefront. In this section, we show that AC is a line, and calculate its slope.

We consider an infinitesimal element, of width Δ*X*, that straddles the wavefront, AC. In figure 4, this element is drawn in a referential system at two different times: *T** and *T***+*Δ*T.* The right side of the element lies in a ‘constant state region’, with *V*_{R}*=**ψ*_{R}*=*1/*T*_{o}*=V*_{o}; the left side lies in an ‘undisturbed region’, with *V*_{L}*=ψ*_{L}*=e=*0. Referring to equation (1.6):Since *ψ*_{R}*=*1/*T*_{o}, *e*_{R}=*e*_{initial}, this can be writtenWe approximate the dimensionless momentum equation withThe velocity of point *X*_{L} increases from 0 to *V*_{o} over Δ*T*. We approximate the acceleration ∂^{2}*U*/∂*T*^{2} as*X*_{L} shifts from the borderline to the constant state region in the time needed to reach the velocity *V*_{o}; the borderline (wavefront) shifts to the left a distance Δ*X* during time Δ*T*. The dimensionless wave speed, *S*, towards the left is thus(3.10)Since *S* is constant along the borderline, AC is a line with slope of 1/*S.*

#### (iii) Distinct regions in the space–time representation

The representative space–time graph of the solution (figure 3) is described by ‘region’ in the following. Note that different choices of parameters might lead to additional classes of regions exhibiting different behaviour, especially in the first few oscillations.

##### Undisturbed regions

As before, triangles like OAC and BFE_{2} contain no wave energy. The proof that *e*=*V*=0 in OAC parallels that for the linear case. The argument is slightly more involved for triangle BFE_{2}, since finding the vertex F requires some care. F is chosen so that all C+ characteristics intersecting the line *X*=1 beneath F emanate from OC or OA and all C+ characteristics intersecting the line *X*=1 above F emanate from above OC. From the initial and boundary conditions on OC and OA, J+=0 on BF. As in the linear case, all C− curves in the undisturbed triangle BFE_{2} emanate from BF. Hence, J−=*ψ*=*e*=*V*=0 in this region. Because *e*=0 in the undisturbed regions, *c*(*e*)=1 is constant and the C+ and C− curves are straight lines with slopes of ±1.

The lower limit of BFE_{2} is the C− line BE_{2}, with a slope of −1. Note that BE_{2} is not a borderline like AC. No C− curves with slope smaller than 1 intersect from above, and no C− curves with slope larger than −1 intersect from below.

##### Constant state regions

ABE_{1}C is a constant state region, in which strain and velocity are constant. This is the region in which all C− curves intersect AB, meaning J−=2/*T*_{o}, and all C+ curves intersect OA or OC, meaning J+=0. Therefore, *V*=*ψ*=1/*T*_{o}. The upper border, BE_{1}, could be any curve between BE_{2} and a line drawn from B with a slope of −1/*c**.

Region FGE_{3}E_{2} contains a constant state sub-section as well (note that the interaction of the wave with the boundary FG is not drawn in figure 3). E_{3}G is drawn such that all C+ curves in FGE_{3}E_{2} pass through CD. In the constant state region within FGE_{3}E_{2}, all C− curves pass through BF. At any point in this constant state region, *ψ*=1/*T*_{o} and *V*=−1/*T*_{o}, meaning strain and velocity are constant. Line E_{2}F demarcates a borderline analogous to AC.

##### Simple wave regions

In regions such as BE_{2}E_{1} which lie adjacent to constant state regions, strain and velocity drop smoothly from their values in the constant state region (on BE_{1}) to zero in the undisturbed region (on BE_{2}). These regions do not exist in the linear case. Since all C+ curves in BE_{2}E_{1} touch the constant state region ABE_{1}C, J+ is constant in this region; hence, *e* and *V* are constant along C− curves, and C− curves in this region are straight lines. The slopes of BE_{2} and D′G′ are −1 and +1, respectively.

Simple wave regions expand over time, as the undisturbed and constant state regions contract and vanish. Upon reflection of the wave, the new simple wave region lies between E_{3}G and D′G′.

##### Boundary regions

Several types of regions occur near the boundaries *X*=0 and 1. The simplest are regions such as CE_{1}D, which are analogous to the constant state boundary regions in the linear case. As in the linear case, point D is chosen so that all C− curves intersecting CD pass through ABE_{1}C; thus, *V*=0 and *ψ*=2/*T*_{o}.

For *μ*≠0, CE_{1} is a borderline like AC, and not a characteristic: C+ curves emanating from CD intersect the steeper C+ curves in ABE_{1}C. E_{1} is the intersection of BE_{1} and the borderline CE_{1}.

Other regions at the boundaries include DE_{1}E_{2}E_{3}, in which waves accelerate from *V=*0 on DE_{1} to their speed on E_{2}E_{3} (strain decreases accordingly), and a region D′DE_{3}, in which energy from the simple wave region is reflected from the boundary.

### (c) Wave broadening

As can be seen in figure 3, the simple wave regions grow over time at the expense of the shrinking constant state and undisturbed regions if *μ*≠0. The wave energy diffuses throughout the specimen until every point in the specimen has some kinetic and/or strain energy at all times. In this section, we re-examine figure 3 to derive an estimate for the time *T*_{f} required for all undisturbed regions to be eradicated.

As described above, the slopes of lines AC and E_{2}F are −1/*S* and 1/*S*, respectively (refer to equation (3.10)) and the slopes of lines BD′ and D′G′ are −1 and 1, respectively. Although the slopes (or average slopes) of curves E_{2}E_{3} and E_{3}D′ are generally different from that of BE_{2}, we neglect this difference when estimating *T*_{f}. This is a reasonable approximation for two different reasons, depending on the size of *T*_{o}. If *T*_{o} is relatively large, so that the end regions are comparable in size to other regions, then the change in the slope is very small, and we may treat E_{2}E_{3}D′ as a linear extension of BE_{2}. On the other hand, if *T*_{o} is small, the whole boundary region is small, and ignoring the now much larger change in slope will not appreciably shift the position of D′.

With this approximation, we examine the simplified space–time plot shown in figure 5, which contains only two types of lines: borderlines with a slope of ±1/*S*, defining the lower boundaries of the disturbed regions, and characteristics with a slope of ±1, defining the upper boundaries of disturbed regions. The ‘height’ of a disturbed region iswhere *X*_{T} is the total distance the wavefront has travelled. The undisturbed region dies out when Δ*T* exceeds the time required for a lower boundary to make a round trip:Solving for *X*_{T}, and noting that *T*=(1/*S*)*X*_{T}, the time *T*_{f} required for all undisturbed regions to vanish is approximately(3.11)As can be seen in figure 5, the final undisturbed region can disappear at an instant when the lower boundary of a disturbed region reaches the end of the specimen. However, the approximation in equation (3.11) does not necessarily find a time that is an integral multiple of 1/*S*. To be more precise, the largest multiple of 1/*S* less than *T*_{f} is the actual time at which the final undisturbed region vanishes from the idealized space–time domain in figure 5.

### (d) Wave frequency

The frequency of the semi-periodic motion in the nonlinear material studied here is a function not only of material properties, but also of the loading rate. This frequency is approximately constant over time, as will be shown by the results of the simulations. We approximate this frequency from the simplified space–time plot shown in figure 5. From figure 5, it is evident that the period of wave motion is dictated by the lower boundary of the active region, whose slope is 1/*S* (equation (3.10)). The period of the cyclical motion is then 2/*S*, and the dimensionless time frequency is approximately

### (e) Wave energy

Closed-form expressions for the work of loading, kinetic energy and strain (potential) energy in a specimen are derived in appendix B. The total work of loading can be written(3.12)Because no dissipative terms exist in the constitutive law, equation (3.12) represents the sum of the total strain and kinetic energy in the specimen at all times. At any point in time, the dimensionless kinetic energy can be found using(3.13)and the dimensionless strain energy can be found using(3.14)As shown in appendix B, dimensional energy is found by multiplying the dimensionless energy by the ‘characteristic energy’, *Γ*:

## 4. Simulations

We developed a finite-difference programme to validate our approximations derived in §3. The code was written in the MATLAB environment, using double precision variables. The differential equation was solved using a central point finite-difference method (Press *et al*. 1992).

We used two different mesh sizes depending on the time range of the simulation. In one group of simulations, where the focus was on steady-state behaviour, we used a larger time range, *T*=[0 100], and a coarser mesh size: Δ*T=*10^{−4}, Δ*X*=10^{−3}. In the other, where the focus was on the first few wave reflections, we used a smaller time range, *T*=[0 2.5], and a finer mesh size: Δ*T=*5×10^{−6}, Δ*X*=10^{−3}. The solution for *U*(*X*,*T* ) was found for all combinations of *μ*={0.01, 0.02, 0.05, 0.08, 0.2} and *T*_{o}={0.01, 0.02, 0.04, 0.07, 0.1, 0.2}.

Noise inherent in numerical schemes caused small fluctuations in the predictions of *U*(*X*,*T* ), which caused errors when taking the derivatives of *U*(*X*,*T* ) needed for calculating the strain, *e*, or the velocity, *V*. To avoid this, *U*(*X*,*T* ) was smoothed by a simple moving window averaging scheme before derivatives were taken (Press *et al*. 1992).

Simulations were checked against the closed-form expression for wave energy (§3*e*) to ensure accuracy.

## 5. Results

In this section, we present results of finite-difference simulations of the problem. In §5*a*,*b*, we check the accuracy of the finite-difference simulations through comparison to exact closed-form expressions, and qualitatively describe the solutions for *U*(*X*,*T* ). In the remaining parts, we validate the approximate formulae developed in §3 by comparing them against the finite-difference results.

### (a) Accuracy of finite-difference results

Accuracy of the finite-difference simulations was checked by (i) refining the grid size until a converged solution was reached and (ii) comparing the total energy from the simulations to the exact work of loading (equation (3.12)). Total energy, calculated at each time increment in each simulation using equations (3.13) and (3.14), remained constant in all simulations after time *T*_{o}. As shown by figure 6, the magnitude of the total energy predicted by finite-difference simulations (grey circles) was very close to the exact analytical expression for the work of loading (solid line) for all loading ramp times *T*_{o}.

### (b) *U*(*X*,*T*) for linear and nonlinear cases

The displacement field *U*(*X*,*T* ) was computed over the length of a specimen (0≤*X*≤1), from the onset of the application of the linear ramp (*T*=0), until *T*=100.

Insight can be gained by examining *U*(*X*,*T* ) both as a function of position at different times (figure 7), and as a function of time at different positions (figure 8). Figures 7 and 8 both show results for a linear ramp applied, over a dimensionless time *T*_{o}=0.1, to specimens with *μ*=0.05.

The ramp, applied at the specimen's right edge (*X*=1), induces a wave to travel from right to left (a ‘departing wave’); the wave reflects from the wall and travels from left to right (a ‘returning wave’).

Figure 7*a* is the displacement distribution at time *T*=2.4, showing the second departing wave for both the linear and nonlinear cases. The linear wave (dashed line) maintains the slope, shape and width (*X*=0.1) of the linear loading ramp, and makes a round trip with a period of *T*=2.0. Wave energy is non‐zero only in the sloped (disturbed) region of the plot; the upper and lower plateaus correspond to undisturbed regions, in which strain and velocity are zero.

The nonlinear wave (solid line in figure 7*a*) spreads over time: the disturbed region is appreciably wider in the inverted returning wave at *T*=3.2 (figure 7*b*) than in the departing wave at *T*=2.4. The nonlinear wave reaches a steady state (figure 7*c,d*), in which the wave spreads to the entire width of the specimen, and the ‘undisturbed region’ is eliminated. The kink in figure 7*c,d* corresponds to interaction with the boundaries.

A second way of viewing these data is to focus on a point and evaluate *U*(*X*,*T* ) as a function of time. This is done in figure 8*a* for point *X*=0.3 and in figure 8*b* for point *X*=0.7. Again, it is evident that the linear wave (dashed line) maintains its original shape, while the nonlinear wave (solid line) spreads. In both cases, energy is initially confined to a region of width *T*_{o}*=*0.1. In the steady state (*T*>10 in this case), this region spreads so that all points contain energy at all times. In the pictured time interval of 20<*T*<24, d*U*/d*T* is always non‐zero, meaning that the two positions shown always have non‐zero kinetic energy.

In figure 8, the linear case repeats perpetually. In each cycle, *U* is zero until the wave arrives. *U* then increases linearly to 1 following the initial ramp, and remains at 1, until the returning wave ramps *U* back down to 0. The amount of time that the displacement stays at 1 increases linearly toward the point of application of the ramp. Only during the ramp-up and ramp-down of displacement are the strain, stress, kinetic energy and strain energy non‐zero.

As predicted by equation (3.10), the material nonlinearity increases the wave speed. This is clear in figure 8: the wave reaches point *X=*0.3 sooner in the nonlinear case, and the frequency of nonlinear steady‐state response (*T* >10) is higher than the linear frequency.

A final feature of the linear and nonlinear steady‐state wave motion that bears mention is symmetry: graphically, one could invert and shift the steady‐state solutions in figure 8*a* for *X=*0.3 and be left with the solution in figure 8*b* for *X=*0.7. Mathematically, this attribute implies that the solutions have two interesting features. The first is that the time average of in the steady state for any position *X* is . Figure 8 shows that at *X=*0.3, this average is and at *X=*0.7, it is . This means that the solutions can be written in the form(5.1)where the time average of *H* is zero. Note that *U=X* is in fact the particular solution to equation (1.5).

The second interesting feature is that the oscillating part of the response *H*(*X*,*T* ) is antisymmetric about *X=*0.5 except for a shift in time of *T*_{p}/2, where *T*_{p} is the time period of the steady-state response:In the linear case, *T*_{p}*=*2; in the nonlinear case, *T*_{p}*≈*2/*S*.

### (c) Transition to steady state and the estimate of *T*_{f}

The transition to steady state is illustrated in figure 9, which depicts the time decay in the peak wave amplitude (difference between maximum and minimum displacement over one cycle) at the mid-point of a specimen with *μ*=0.01. Circles are peak wave amplitudes, from the simulation results. The transition to steady state occurs more quickly when steeper loading ramps are applied. Also shown in figure 9 (dotted lines) are the estimates of steady‐state transition times, *T*_{f}, from equation (3.11), which proved to be reasonable estimates for all cases tested.

### (d) Space–time map of strain and the estimate of wave speed

A space–time map like that shown in figure 3 is plotted in figure 10, using the results of a finite-difference simulation of a nonlinear case (*T*_{o}=0.2, *μ=*0.08). Figure 10 shows only the strain, *e*, calculated from the *U*(*X*,*T* ).

As in figure 3, the white regions of the plot are undisturbed regions where displacement is constant (*U=*1 in the undisturbed regions with an edge along *X=*1; *U=*0 in the undisturbed regions with an edge along *X=*0). The shaded regions are disturbed regions, where strain is non‐zero; darker shading indicates higher strain. The thickness of this region increases with time, as energy diffuses throughout the specimen.

Three lines in figure 10 are drawn from point (*X*,*T* )=(1, 0), where the wave starts. The lower dashed line, with a slope of −1/*c**, is the same as AC″ in figure 3 (*c** is defined in equation (3.9)). The upper dashed line, with a slope of −1, is the same as line AC′ in figure 3. As described in §3*b*(i), the wavefront can be any line between (and including) the two dashed lines. The dark line between the two dashed lines has a slope of −1/*S*, where *S* is the estimated wave speed (equation (3.10)). In this case, the estimate is very close to the actual borderline.

Two other solid lines are drawn in figure 10: the first, with a slope of −1, emanating from the point (*X*,*T* )=(1,*T*_{o}=0.2), and the second, with a slope of +1, emanating from the first's end point (*X=*0, *T=*1.2). These two lines are BD′ and D′G′ in figure 3. These lines provide a good estimate of the upper border of the disturbed region.

The sub-regions of the space–time plot in figure 3 can all be seen in figure 10. A constant state region is visible in the lower part of each disturbed region. Simple wave regions are visible on ‘top’ of these constant state regions (that is, at a greater time *T* for a particular position *X*, or a value of *X* further from the wavefront for a particular time *T*). Boundary regions are visible near *X=*0 and *X=*1. The triangular region CDE_{1} (refer to figure 3) is evident, where *ψ*(*e*)=2/*T*_{o} and *V=*0. Other boundary regions are visible above this triangle, as in figure 3.

### (e) Wave speed

The closed‐form approximation for wave speed (equation (3.10)) was evaluated as a function of *μ* and the loading time, *T*_{o}. Although waves in nonlinear media spread over time until they reach a steady state, the simulations support the approximation that the time period remains almost constant over all wave cycles, which was made when deriving equation (3.10).

Figure 11 compares the predicted first cycle wave speed (which, again, is very close to the steady‐state wave speed), to the estimates of equation (3.10), for several values of *μ* and *T*_{o}. The grey circles, representing the wave speed from the simulations, match the black lines, representing predictions using equation (3.10), to within about 5%.

### (f ) Energy partitioning

A central result of this work is the ability to determine loading rates at which dynamic effects become important. We evaluate this by considering the fraction of the energy input to the specimens that takes the form of kinetic energy; the lower the fraction of kinetic energy, the smaller the importance of dynamic effects. Figure 12 compares the kinetic energy fraction in a nonlinear specimen, calculated from equation (3.13), to that in a linear material, which can be calculated in closed form (appendix C). In both cases, ‘nodes’ appear for ideal loading rates, in which dynamic effects disappear entirely. These nodes correspond to loading at a rate such that the completion of the loading ramp coincides with the return of the wavefront to the loading point. Results for lower values of *μ* fall between the two curves pictured. In all cases, loading ramps applied over *T*_{o}>2 result in dynamic effects of less than 5%.

## 6. Discussion

A rapid stretch applied to one end of a relaxation specimen instigates a cyclical wave. In linear elastic media, this wave sustains its shape; in nonlinear media, this wave eventually diffuses over the entire length of the specimen, until it reaches a steady state.

We derived expressions for the time required to reach this steady state, for the wave speed and for the work of loading. Comparison with finite-difference simulations showed these to be accurate in the parameter ranges of interest. The wave speed estimate was accurate to within 5%, with the greatest deviation for the shortest ramp times, *T*_{o}. The source of this discrepancy is the approximation made in deriving equation (3.10). The approximation that the velocity of a point on the wavefront increases linearly from 0 to *V*_{o}(∂^{2}*U*/∂*T*^{2}=*V*_{o}/Δ*T* ) is acceptable only for small or moderate nonlinearity and longer ramp times.

We conclude with a discussion of our decision to model a step stretch with a quick ramp and a discussion of how to apply the results presented in this paper to design and interpret relaxation tests on nonlinear media.

### (a) Limitations as the loading ramp approaches a step

We modelled an ideal ‘step’ stretch with a linear loading ramp. As discussed in this section, this was necessary, because no solution exists for a ‘step’ loading of nonlinear elastic materials of the kind represented by equation (1.3). When an ideal step stretch is applied to one end of the specimen, *U* is discontinuous at *X=*1 and *T=*0, and its partial derivatives with respect to both *T* and *X* directions approach infinity as a Dirac delta function and . Take second derivatives(6.1)and(6.2)where *δ*′(*x*) is defined throughFrom this definition, *δ*′(*x*) must be zero everywhere except *x=*0, where it approaches both positive and negative infinity.

Substituting equations (6.1) and (6.2) into equation (1.5) at *X=*1 and *T=*0:(6.3)For the linear case, *μ=*0, yieldingBoth sides of this equation are infinite, but the equality is valid if *A*_{T}=*A*_{X}.

However, for *μ*≠0, the right side of equation (6.3) is one order of *δ* greater than the left side. The equality cannot be satisfied in this case, meaning that no solution exists for an ideal step stretch applied to a nonlinear material following equation (1.3).

### (b) Design of quasi-static relaxation tests

Force measurements from a quasi-static relaxation test can be interpreted to characterize a specimen's constitutive behaviour. If the strain field in such a test is uniform and constant throughout a specimen, the relationship between stress and strain shows the constitutive law. Figure 12 suggests two approaches in achieving quasi-static loading conditions.

The first approach is to ‘tune’ the loading time, *T*_{o}, to eliminate all post-loading wave motion: for several values of *T*_{o}, no kinetic energy remains after the first wave cycle. This is the case when the curves in figure 12 intersect the *T*-axis. In a linear specimen, this occurs when *T*_{o} is an integral multiple of 2; in a nonlinear specimen, this occurs for slightly lower values of *T*_{o}. An additional benefit of tuning the loading time, *T*_{o}, is that, if *τ* is known and if *T*_{o} is tuned to the smallest value that eliminates kinetic energy, then *μ* may be estimated by solving equation (3.10):This approach to estimating *μ* is extremely sensitive to measurement error, however. A 1% error in the wave speed measurement converts to a 25% error in *μ*.

The second approach is to choose *T*_{o} sufficiently large that an acceptably small fraction of the work of loading becomes kinetic energy. Figure 12 shows that as the loading time increases, the wave disturbance decreases on average. The ratio of kinetic energy to total energy will be smaller than 2% for *T*_{o}>4, or *t*_{o}>4*τ*.

#### (i) Case study: pure collagen

As an example, we assess the pure collagen specimens tested by Pryse *et al*. (2003). In these specimens, the time-scale, *τ*, was about 5 ms, and the nonlinearity coefficient from equation (1.3) was *ca* *k=*1.2. The length of their specimens was *l*_{o}=15 mm, and the displacement amplitude was *u*_{o}=0.3 mm. These combine to yield *μ*=0.05. Referring to figure 12, dynamic effects will be smaller than 5% for *T*_{o}>2 or *t*_{o}>2*τ*. Since 5% is an acceptable error in this modelling process, a test could be considered quasi-static for loading speeds *u*_{o}/*t*_{o}<0.03 m s^{−1}.

### (c) Interpretation of dynamic tests

Fitting the constitutive law, equation (1.3), to a material involves identifying *E*_{o} and *k*. We present two approaches for deriving these from results of dynamic experiments. First, *E*_{o} and *k* can be found by (i) measuring the cycle time, *T*_{p}*≈*2/*S*, in two tests with different ramp times, *t*_{o}; (ii) fitting the results to the expression for *S*; then (iii) solving for *τ* and *μ*. This requires very fast loading of the specimen for acceptable accuracy.

The second approach involves analysing the stress and strain averaged over one period, and . and decay down to steady-state values. The goal is to understand how the relation between and can explain the constitutive law.

Equation (5.1) showed that the steady‐state displacement of any point in the specimen can be written as a summation of its position, *X*, and a periodic function whose time average is zero, *H*(*X*,*T* ). The dimensionless strain at any point in the specimen iswhere ( )′=∂/∂*X*. Since strain must exist everywhere, *H* must be continuous and differentiable with respect to *X*, and *H*′(*X*,*T* ) must be a periodic, piecewise, continuous function whose average over each time period is zero. Hence, the steady‐state time average of *e*(*X*,*T* ) is equal to one everywhere in the specimen.

Stress can be found using equation (1.6):The time average of stress over one period is(6.4)where *T*_{p} is the period (*T*_{p}≈2/*S*). The average stress depends on the mean square of the oscillatory part of the strain, which can be shown to be independent of *X*.

It is always possible to choose *T*_{o} sufficiently large that the integral in equation (6.4) is negligibly small, since the mean square of *H*′ scales with the kinetic energy. Therefore, in ‘quasi-static’ conditions (refer to §6*b*), equation (6.4) reduces toor in dimensional form(6.5)*E*_{o} and *k* may thus be found from the average stress measured in two or more tests, each with a different stretch *u*_{o}.

### (d) Viscoelastic materials

We assumed our specimens to be elastic. However, relaxation testing is used mainly for viscoelastic materials. The first approach to interpreting dynamic experiments presented in §6*c* is unaffected by moderate viscous damping, and can be used to estimate *E*_{o} and *k* provided that waves persist for one complete cycle.

If the wave period is much shorter than the relaxation time constant, so that viscoelastic dissipation is small in the first few wave cycles, then approximating the material as initially nonlinear elastic is reasonable. The second approach for estimating *E*_{o} and *k* is valid in such circumstances. However, after a time period on the order of the shortest characteristic ‘relaxation time’ of a material (Flügge 1975), the effects of viscous dissipation must be addressed by an extension of the work in this paper. Techniques from this are firmly in place for the semi-infinite case (Moodie & Barclay 1991; Oncu & Moodie 1993).

## Acknowledgements

We are grateful to the National Institutes of Health for providing financial support for this work under grants AR47591 and GM038838. We also thank Niloufar Ghoreishi, Juan Pablo Marquez and Tony Pryse for their helpful comments and suggestions.

## Footnotes

- Received July 22, 2004.
- Accepted November 2, 2004.

- © 2005 The Royal Society