Modeling of shock propagation and attenuation in viscoelastic components

Protection from the potentially damaging effects of shock loading is a common design requirement for diverse mechanical structures ranging from shock accelerometers to spacecraft. High damping viscoelastic materials are employed in the design of geometrically complex, impact-absorbent components. Since shock transients are characterized by a broad frequency spectrum, it is imperative to properly model frequency dependence of material behavior over a wide frequency range. The Anelastic Displacement Fields (ADF) method is employed herein to model frequency-dependence within a time-domain finite element framework. Axisymmetric, ADF finite elements are developed and then used to model shock propagation and absorption through viscoelastic structures. The model predictions are verified against longitudinal wave propagation experimental data and theory.


Introduction and background
Mechanical parts are often required to withstand shock loads.Viscoelastic materials exhibit high damping characteristics and are commonly used in the design of impact-absorbent components.Shock transients have a wide frequency spectrum, thus it is imper- * Corresponding author: Dr. Razvan Rusovici, STI Technologies Inc., 1800 Brighton-Henrietta Townline Rd., Rochester, NY 14623, USA.Tel.: +1 716 424 2010; Fax: +1 716 272 7201; E-mail: rrusovici@sti-tech.com.ative to properly model the frequency dependence of material behavior, such as the storage and loss moduli.
The classical theory of elasticity states that for sufficiently small strains, the stress in an elastic solid is proportional to the instantaneous strain and is independent of the strain rate.In a viscous fluid, according to the theory of hydrodynamics, the stress is proportional to the instantaneous strain rate and is independent of the strain.Viscoelastic materials exhibit both solid and fluid characteristics [1,2].Such materials include plastics, rubbers, glasses, ceramics, and biomaterials.Viscoelastic materials are characterized by constant-stress creep and constant-strain relaxation.Their deformation response is determined by both current and past stress states or, conversely, the current stress-state is determined by both current and past deformation states.It may be said that viscoelastic materials have "memory"; this characteristic constitutes one foundation on which their mathematical modeling may be based [1].In polymers, material damping is a direct result of the relaxation and recovery of the long molecular chains after stress.

Linear viscoelasticity
Linear viscoelastic materials may be defined by either differential or integral constitutive equations [1,2].Only the differential constitutive equations are reviewed here.The differential form of the constitutive law for a one dimensional linear viscoelastic solid is shown in Eq. (1).
The quantities σ and ε are the one-dimensional stress and strain, respectively; P and Q are the differential operators, which are further expressed as (2) where p j and qj are constants.The number of time derivatives retained in the operators P and Q are denoted by m and n, respectively; j is the order of differentiation.
The Maxwell model of an one-dimensional viscoelastic material consists of a linear spring and a linear dashpot connected in series.The Kelvin model consists of a linear spring and a linear damper connected in parallel.To better approximate material behavior over a frequency range, more than one spring or linear damper may be included in the model.The three-parameter model of a standard viscoelastic solid consists of a linear spring in series with a linear Kelvin element, or a Maxwell element in parallel with a mechanical spring For harmonic forcing and response, the differential linear viscoelastic constitutive law, Eq. ( 1), yields where K * is the complex material modulus; K and K are the material storage and loss moduli, respectively, corresponding to the circular frequency ω.Equation (3b) states that, at a given frequency, a phase shift exists between an oscillatory stress and its corresponding displacement response.The effective moduli of linear viscoelastic materials are functions of frequency and temperature [3,4]; however, they are independent of stress and strain.The frequency dependence of the complex modulus components must be captured accurately; simple models, such as Maxwell or Kelvin, are unable to do that.The Maxwell model behaves like a fluid at low frequencies, while the Kelvin model becomes infinitely stiff at high frequencies.

Internal variable modeling
The complex modulus variation with frequency may be approximated using the concept of internal variables.Various modeling techniques may be used to incorporate linear, frequency-dependent material behavior in structural dynamics, and, in particular, Finite Element (FE), analysis.A few of these methods are mentioned here.
McTavish and Hughes extended the Golla-Hughes model and formulated the GHM (Golla-Hughes-McTavish) model for linear viscoelastic structures [6][7][8].In this formulation, the material modulus in the Laplace domain is modeled as the sum of minioscillators, consisting of a viscous damper, spring and mass.The motion of the mass represents the internal dissipation coordinate.This method leads to a system of second order differential equations of motion in the frequency domain, where the mass, stiffness, and damping matrices are augmented by the internal dissipation coordinates.A difficulty encountered with this particular method is that the stiffness matrix is ill conditioned and spectral decomposition is needed to remedy the difficulty.
Yiu [9,10] introduced another first-order technique to develop FE that includes frequency-dependent damping behavior of linear viscoelastic materials.This technique employs a generalized Maxwell model to represent material behavior.The dynamic stiffness matrix is obtained in the Laplace domain.Internal coordinates are needed to convert frequency domain equations to the time domain.Each internal coordinate is related to a Maxwell element.Later, the viscoelastic operator was modified by the inclusion of a dashpot in parallel with the generalized Maxwell model, in order to rule out instant elastic response.As with the GHM method, the Yiu method only models linear viscoelastic materials.
Lesieutre [11,12] and his co-workers presented the "Augmenting Thermodynamic Fields" (ATF) approach to model frequency-dependent material damping of linear, viscoelastic structures in a finite element context.The ATF models the dissipative behavior of linear damping materials.The ATF evolution equation is determined from the irreversible thermodynamics assumption that the rate of change is proportional to its deviation from an equilibrium value.This results in coupled partial differential equations in terms of the displacement field and the gradient of the ATF.The material constitutive equations are derived from the Helmholtz free energy function.A first-order finite element model is obtained from the discretized system of equations.The coupled equations of motion are discretized using the method of weighted residuals, and a first order equation system is obtained.The ATF parameters are found by iteratively curve-fitting the predicted complex modulus to experimental data.Lesieutre [13] developed the "Anelastic Displacement Fields" approach as an extension of the ATF method.This method considers the effect of material anelasticity on the displacement field, as opposed to directly modeling physical damping mechanisms.Lesieutre and Bianchini [14] developed rod, beam and plate, ADFbased finite elements.Govindswamy [15] used ADF rod elements to study longitudinal wave propagation along a viscoelastic bar.
In this work, the ADF method is employed to develop axisymmetric and plane stress finite elements that are capable of modeling the frequency dependent behavior of linear viscoelastic materials.The newly developed finite elements may be used to model and analyze behavior of structures subjected to shock loads, such as propagation of longitudinal waves through viscoelastic rods, and shock absorbent mechanical filters.

Longitudinal wave propagation through rods
The wave equation governing the one-dimensional longitudinal rod motion is where E is the Young's modulus of elasticity, which may be frequency dependent, ρ is the material density, u(x) is the displacement at a station x along the rod and q(x) is a body force acting along the bar [16,17].The hypothesis behind the equation is that during deformation, initially plane sections of the rod remain plane.This equation assumes uniaxial strain and neglects the effects of the Poisson's ratio on the lateral displacements.In an elastic rod the Young's modulus is frequency-independent, and the propagation velocity is Propagation or bar velocities in most metals are around 5000 m/s.The particle velocity, in an elastic bar, is defined by [16] V where σ is the stress in the bar at the given time and location, ε is the corresponding strain and c 0 is the bar velocity.In contrast to propagation velocities, particle velocities in bars are several orders of magnitude less, since the stress σ is typically several orders of magnitude smaller than the modulus of elasticity E. Consider that an incident stress pulse reaches the fixed boundary of a fixed-free rod.At that boundary, the reflected stress pulse has the same sign as the incident stress pulse.At the free boundary, the pulse changes sign.On a freefree rod, the initial stress pulse reflects repeatedly due to end reflections.Velocity doubling occurs at the free ends of the bar.All points on the bar exhibit "a series of jerky movements" [16], as seen in Fig. 1.
During longitudinal wave propagation through a bar, a phenomenon called geometrical dispersion is encountered.As a result of the longitudinal displacement and Poisson's ratio ν, displacements in the other two coordinates y and z occur (v and w, respectively): The wave propagation equation, which accounts for lateral inertia effects, was developed by Love [17].Geometrical dispersion becomes a factor only when the wavelength of a given frequency component present in the incident stress pulse Λ is on the same order as, or smaller than, the bar radius r Lateral inertia causes high frequency components of the pulse to travel at lower velocities than lower frequency components [16].When a longitudinal pulse travels through a viscoelastic rod, both viscoelastic attenuation and dispersion, and geometric dispersion phenomena may affect it.Since viscoelastic materials exhibit frequency dependence of their material moduli, the propagation velocity also varies with frequency [18,19].When a longitudinal stress-pulse propagates along a viscoelastic bar, high-frequency components travel faster and attenuate more rapidly than low-frequency components; note that this effect is opposite to that produced by geometric dispersion.

Experimental setup
The setup used to investigate longitudinal stress pulse propagation in a PERSPEX TM viscoelastic rod is described next.The bar is simply supported, but not constrained axially at two longitudinal locations.A compressed air-driven DELRIN TM bullet, hits the bar at one of the free ends.The longitudinal impact creates a stress wave that propagates along the bar.The compressive stress pulse generated by the projectile impact usually resembles a half-sine or versed cosine shape.Velocity, shape (blunt or sharp impact end), and material of the projectile are all factors that influence the duration and shape of the resulting stress pulse.The setup is similar to the one used for shock accelerometer calibration with a Hopkinson bar [20][21][22], and is shown in Fig. 2.
The impact force is measured with a force gauge and the strain in the middle of the bar is measured with a pair of strain gauges.Two strain gauges are needed to compensate for any bending strains.A Gage TM data acquisition board with 1 MHz sampling rate, records impact force and strain gauge data.The force and strain gauge data is recorded and processed by Labview code adapted from the Hopkinson bar calibration code developed in part by the author [22].Then, the force data is used as input to the finite element code, and the strain gauge record is compared to the model prediction.

Finite element development
This section shows how frequency-dependent material behavior modeled using the Anelastic Displacement Fields (ADF) is incorporated in a finite element framework, and presents the development of the triangular axisymmetric finite element [23].First, the linear ADF equations of motion and the material constitutive equations are presented; then, the ADF finite element equations are developed.
According to the theory of elasticity for the axisymmetrical problems, two components of the displacement completely define the stress state: u, the displacement in the radial direction r, and v, the displacement in the longitudinal direction z.A radial displacement induces a strain in the circumferential direction.Thus, the following are true for axisymmetric problems [24] σ r , σ θ , σ z , τ rz = 0 ε r , ε θ , ε z , γ rz = 0 (9) The axisymmetric hypotheses also apply to the anelastic stresses and strains:

ADF modeling
The structural response of a system is determined in computational structural dynamics by solving a set of n simultaneous, time-dependent system of equations where M, K and C are the mass, stiffness and damping matrices, respectively and {f (t)} is the time dependent applied load.Lesieutre [12] introduced the ADF method to model damping in a finite element framework.This method considers the effect of material anelasticity on the displacement field: the displacement field u consists of an elastic component and an anelastic components: u E (x, t) and u A (x, t) respectively = u E (x, t) where the anelastic component may be written as the sum of n components, as described by Eq. ( 12).In the ADF model, stress-strain relationships are written in tensor form as [13,14] where ε A are the anelastic strains.The anelastic stresses are Finally, the coupling between the elastic and anelastic field is where Ω is the inverse of the relaxation time, and E A ijkl are anelastic material constants.
The equations of motion are where f i are the body forces, and ρ is the density.More than one ADF field may be used to for modeling purposes.The frequency dependence of an individual scalar modulus, β * m , is characterized by [13] where ω is the circular frequency (rad/s), β 0 m is the complex modulus, β mr is the relaxed or low-frequency modulus, Ω n is the characteristic relaxation time at constant strain corresponding to the ith anelastic displacement field, ∆ mi is the relaxation strength corresponding to the ith ADF field and modulus m.The high frequency, unrelaxed, modulus is [13] and The constitutive parameters C i show the coupling between the total displacement field and the physical relaxation process The ADF model input parameters are determined by curve fitting the predicted material moduli to corresponding experimental values.The input parameter values are adjusted iteratively to minimize differences between predicted and experimental values in a least square sense.The number of anelastic displacement fields used usually determines how well the predicted material data compares to the actual value over a given frequency range.

Finite element development
Simple triangular axisymmetric ADF finite elements are developed.Each node has two degrees of freedom (DOF): u (the displacement in the radial direction r) and v (the displacement in the vertical direction z).The element area is denoted by ∆A.The displacements are approximated using the following shape functions at each node (denoted respectively by subscripts e, f and g) of the triangular finite element By substituting the shape functions, the strain matrix is obtained as a function of the displacements δ and the anelastic displacement vector δ A : and The strain energy is After integrating over the volume, the element stiffness matrix becomes

) [B]rdrdz
The element damping-matrix is found using similar considerations [13] [c] = [B]rdrdz The mass matrix of the ADF axisymmetric element does not include any contribution from the anelastic displacement fields and resembles the mass matrix of a simple, elastic, triangular axisymmetric finite element Thus, the stiffness, damping and mass matrices are developed for a simple, triangular, ADF axisymmetric element.

Results
In this section, model predictions are compared to longitudinal wave propagation experimental results.Two cases are analyzed: the longitudinal wave propagation along a viscoelastic rod and the modeling of mechanical filters subjected to shock.

Wave propagation through a viscoelastic rod
A shock force usually excites frequency components spread over several decades [23].Thus, it becomes necessary for analysis purposes to define a high cutoff radial frequency, which also allows for reasonable solution accuracy [25,26].The goal is to determine that cutoff radial frequency ω c0 and to define a corre- sponding finite element mesh.The cutoff frequency determines the critical wavelength L w .
An axisymmetric, three-ADF finite element model is employed to simulate the transient response of a 1.8-meter long, free-free PERSPEX bar subjected to a half-sine, longitudinal impact force.The bar has a 15.7-mm diameter.The mesh consists of 250 nodes in the longitudinal direction by two nodes in the radial direction.The input force has a duration τ of 150 µs.The Fourier amplitude spectrum of a half-sine pulse becomes null at a frequency f , where In the experiment, the impact force lasts approximately 100 µs.This corresponds to a cutoff frequency of 15 kHz.The total time for the corresponding wave to travel a distance equal to an element length is where c is the wave velocity at the corresponding frequency (c for PERSPEX TM at 15 kHz is 2193 m/s).
Assuming that n1 time steps are necessary to properly represent the waveform of wavelength L w , the time step is where n1 is recommended to be equal or greater than 20 in finite element practice [25,26].The element length is The time step chosen is 1 µs.The time step is chosen to be smaller than the calculated values to allow for better solution accuracy.The cutoff frequency also determines the frequency range of interest for ADF modeling of material parameters.The Newmark method is used here to numerically integrate the equations of motion, since this particular method is inherently stable.
The parameters presented in Table 1 were obtained through a least square curve-fit of the PERSPEX TM Young's modulus and loss factor data using a three-ADF model.In practice, material modulus-versusfrequency data often varies from one batch to the next.The ADF parameters can be slightly adjusted to better match displacement time response.
The relaxation strengths of the shear ∆ G modulus and bulk modulus ∆ K are assumed to be equal.The Young's modulus is approximated better than the loss factor over the frequency range of interest.More ADF would increase accuracy, but would require additional computational resources.Figure 5 shows the three-ADF model prediction versus the experimental strain time record in the middle of the free-free PERSPEX TM bar.The time record starts at the beginning of the first pulse, and only the first two pulses are considered.Figure 5 also shows the viscoelastic attenuation of the initial pulse, as it travels along the bar, a phenomenon captured well by the finite element model.The model predicts the peak value of the first pulse to be within 9% relative error to the experiment, while the second pulse is predicted within 4% relative error.The frequency content of the predicted second pulse shows slightly less dispersion than in reality.That is because the three-ADF material parameter approximation is valid over a limited frequency range.Note in Fig. 6 the decay of the mechanical energy in the viscoelastic bar.The mechanical energy is the sum of the kinetic and strain energies.The strain energy component decays and the total mechanical energy asymptotically approaches the constant-velocity kinetic energy of the bar.Such results could not be obtained using a purely elastic finite element model.The preceding experiment was repeated for a Perspex bar with a 28-mm diameter and 1.9-m length.The predicted strain in the middle of the bar is compared against the experiment in Fig. 7.
Again, the model predicts the peak value of the first pulse to be within 10% relative error to the experiment, while the second pulse is predicted within 3% relative error.sion.Geometric dispersion occurs when the diameter of the bar is on the same order of magnitude as, or longer than the wavelength of frequency components contained in the stress pulse.The stress pulse is dispersed as a result of the effect of lateral inertia.Geometric dispersion causes high frequency pulse components to travel slower than lower frequency pulse components [16].As a result, high-frequency oscillations follow the displacement step pulse associated with longitudinal wave propagation along bars.These high-frequency oscillations that follow a step displacement pulse, indicative of geometric dispersion, were recorded and published by Kolsky [19].

Modeling of geometrical dispersion in a PERSPEX
The next numerical simulations show that geometrical dispersion is predicted by an ADF based finite element model.In these simulations, a free-free PERSPEX TM bar is hit at one end by a half-sine force of 1368 N magnitude and 4.35 kHz central frequency.In one simulation the bar has a 15-cm diameter, and in the other simulation the bar has a 1.5-cm diameter.The bar is 1.8 m long.The FE mesh that models the bar has 200 nodes in the longitudinal direction by 2 nodes in the radial direction; the element length is 0.75 cm.A single-ADF, triangular, axisymmetric FE is used.The time step chosen is 2.5 microseconds.This numerical arrangement allows good modeling of frequency components of up to 14.5 kHz in the response.Pulse dispersion is expected in the 15-cm diameter bar, since the wavelength corresponding to a 14.5 kHz frequency is on the same order of magnitude as the larger diameter.
Figure 8 compares the predicted bar end displace-ments for the two simulations.Each end displacement is normalized to its maximum value.The two bars have different masses, but the same force acts upon them: as a result, the magnitudes of the corresponding end displacements are different.Note in Fig. 8 that a series of high-frequency oscillations follow each sharp rise of the end displacement corresponding to the 15-cm diameter bar, which are the effects of geometrical dispersion.The geometrical dispersion effects are not observable in the end displacement of the smaller diameter bar, however.The effects of viscoelastic attenuation are present in both displacement traces and are represented by slope changes in the displacement time record. Figure 9 shows the mechanical energy in the 15-cm bar, which is the sum of the kinetic and strain energies.The strain energy component decays and the mechanical energy asymptotically reaches the constant-velocity kinetic energy of the bar.

Modeling of mechanical filters subjected to shock
Mechanical filters are often employed in shock accelerometers to protect the sensing element from shock transients [27][28][29][30].It is of interest to quickly determine whether a given material will provide the desired damping characteristics.
To verify this, a mechanical filter made of Buna N B-1 [31,32] rubber was glued to the end of a 2.03-m long Hopkinson bar.The filter was a 3-mm thick disk with a 1.9-cm diameter (same diameter as the bar).A 30-gram mass was glued on to the filter; the mass had a shock accelerometer screwed into it.The shock accelerometer did not have a mechanical filter itself.A force gauge was screwed into the opposite end of the bar.A DELRIN TM bullet driven by compressed air hit the force gauge.The impact force (shown in Fig. 10) and acceleration were recorded using the experimental setup described in Section 2 and Fig. 1.The force data was used as input into the finite element code; the acceleration record was integrated with respect to time and the resulting velocity record was compared to the model prediction.The ambient temperature during the experiment was approximately 25 • C.
Axisymmetric, single-ADF finite elements were employed to model the filter, while elastic axisymmetric finite elements were used for the bar.The bar was modeled using a mesh made of 300 nodes longitudinally by 2 nodes radially.The rubber disk was represented using a mesh formed by 10 nodes longitudinally by 2 nodes radially.To allow proper representation of frequency components up to 6 kHz, an integration time step of 1 microsecond was chosen.
The following formula to compute the relaxation magnitude for a single-ADF finite element [14]: where η p is the peak loss factor.The inverse of the relaxation time is where ω p is the circular frequency corresponding to the peak loss factor value.All the ADF parameters used are shown in Table 2.
The predicted longitudinal particle velocity and experimental velocity of the mass versus time are compared in Fig. 11.The magnitude of the first velocity peak is predicted within 5% relative error to the experimental value.
Figure 12 shows the predicted longitudinal displacement of the mass over time.Note that filter resonances are superimposed on the displacement steps.Figure 13 depicts the mechanical energy in the bar versus time.
The filter absorbs the energy as the pulse passes through it and causes the steps observed in the energy time record.A new analysis was performed using three-ADF finite elements to model the filter, since it was shown in Reference [13] that using more ADF coordinates improves the accuracy of the material model.The ADF parameters were found this time by curve-fitting the predicted moduli to measured data found in Reference [32].Figures 14(a) and (b) show a comparison of the three-ADF curve-fit of the storage and loss modulus, respectively, to corresponding experimental data.
The longitudinal particle velocity predicted using a three-ADF (see Table 3) axisymmetric finite element model is compared against the experiment, as shown in Fig. 15.
The magnitude of the first velocity peak is also predicted within 5% relative error to the experimental value.

Conclusions
In conclusion, the experiments and simulations summarized here verify that ADF axisymmetric finite elements provide a good model for capturing the behavior of viscoelastic components under shock loading.Such models are capable of capturing wave propagation phenomena, including geometric dispersion and viscoelastic attenuation of longitudinal waves in rod, and behav-ior of mechanical filters in realistic shock conditions.ADF three-dimensional finite element models may thus be successfully employed to design mechanical filters for various shock-isolation applications.

Fig. 1 .
Fig. 1.End and middle displacement on free-free bar.
Figures 3 and 4 show the three-ADF model curve fit of the loss modulus and Young's modulus, respectively, versus the corresponding experimental values.

Table 1
ADF parameters used for PERSPEX TM

Table 2
Single-ADF parameters used for Buna N B-1 rubber ∆ G = ∆ K ; the relaxation strengths corresponding to the shear and bulk modulus, respectively, are considered to be equal. *

Table 3
Three-ADF parameters used for Buna N B-1 rubber