Printer Friendly

A finite element simulation to longitudinal impact waves in elastic rods.


Investigation of wave propagation in a rod due to impact has a long history. Bernolli, Navier, Poisson, and St. Venant are among the great researchers who investigated this problem. Good reviews of the treatments of longitudinal waves in rods produced by impact are offered in [1, 2]. Recently, due to the presence of powerful computers, new computational methods are applied to solve this classical problem. An analytical simulation, symbolic solution, and a solution using time delay method are developed [3, 4, 5]. Both theoretical and experimental researches were conducted [6, 7, 8] and a review of the experimental studies is offered in [9].

One practical devices utilizes longitudinal wave propagation in rods is Hopkinson (or Davies) bar. The Hopkinson pressure bar is a long, thin, and elastic rod, in which a stress wave is generated at one end by a projectile impact. The projectile is a rigid mass or a striker bar. At the other end of the bar the generated wave can be used in many applications. The propagation of shock wave in Hopkinson bar, used to calibrate shock accelerometers under high acceleration levels and a wide frequency bandwidth, is modeled [10, 11]. Insertion of a deformable disk between the projectile and the bar can decrease the wave dispersion, hence, a commercial finite element code is utilized to investigate dispersion in the bar and to find the optimum characteristics of the inserted deformable disk [12].

Wave propagation can be used in determination of mechanical properties of materials. Some dynamic strength material constants are obtained using the split Hopkinson pressure bar [13]. The split tensile Hopkinson bar tests are interrupted to evaluate the damage in the materials at high strain rate [14]. The evaluation of the coefficient of restitution, through numerical simulation of impact of a rigid mass and a slender elastic rod, is investigated [15, 16]. Furthermore, there is an increasing interest in using wave propagation in crack detection, for example, wave propagations in cracked beams and plates are examined [17,18].

Some machine elements are rod-like bodies that are subjected to impact loading during their functional operations. Examples are encountered in piling, percussive drilling and hydraulic hammering. Due to the elasticity of these axial elements, waves propagate through them while they are in translational motion. In the same time, it is obvious that wave propagation is gaining more potential in determination of mechanical properties of materials and non-destructive testing methods. Therefore, reliable finite element models are needed to be used in the simulation of the propagation of waves. Though, some of the previously reviewed works used finite element models in their analysis, the impact forces were assumed, see [17, 18], or calculated using methods that are highly time consuming, see [15, 16]. In this paper, the contact force is calculated using an efficient approach. A finite element model is constructed to represent impact of a rigid mass on a flexible rod. The model utilized St.Venant classical impact model. The two cases of free-free and free-fixed elastic rods are investigated. A numerical scheme is formulated depending upon Newmark implicit time stepping method and Newton-Raphson iterative method. The contact force is calculated and the wave propagation in the rod is simulated. To enhance the understanding of the complicated physical phenomenon, a simulated visualization of the propagation of the impact wave through the bar is monitored.

Mathematical Modeling

It is assumed that the rod has massm, Young's modulus E, density [rho], cross-sectional area A and length l. The rod is initially at rest and is struck on the right end x = l at the initial time t = o by a moving rigid mass [m.sub.0] with initial velocity [v.sub.0]. The displacement of the rigid mass at time t is donated by q(t) and the displacement of the rod at position x and time t is given by u(x, t), as in Fig.1.


The governing equation for the longitudinal wave in the rod is

[partial derivative]u(x, t)/[partial derivative][t.sup.2]=[c.sup.2][[partial derivative].sup.2]u(x, t)/[partial derivative][x.sup.2] (1)

where c is the wave propagation velocity

c = [square root of E/[rho]] (2)

The strain [epsilon](x,t) in the rod is given by

[epsilon](x,t) = [partial derivative]u(x, t)/[partial derivative]x (3)

For elastic rod, the stress is proportional to strain, or

[sigma](x, t) = E [partial derivative]u(x, t)/[partial derivative]x (4)

As the contact established between the mass and the rod, both the mass and the contact end of the rod (x = l) are assumed to have the same velocity [v.sub.0]. Therefore a compression wave is created in the rod. The wave travels along the rod and reflected at the other end (x = 0). During contact period, displacement q(t) and velocity [??](t) of the mass are the same as those of the contact end of the rod (x = l).

q(t) = u(l, t) and [??](t) = [partial derivative]u/[partial derivative]t(l, t), 0 < t < [t.sub.c] (5)

where [t.sub.c] is the contact period.

The contact persists as long as the contact force between the mass and contact end of the rod does not vanish. The contact force equals the stress at the contact end times the rod's cross-sectional area i.e.

F(t) = EA [partial derivative]u(l, t)/[partial derivative]x (6)

The motion of the rigid mass is governed by

m d[??]/dt = F(t) (7)

Equations (1), (6), and (7) are the equations of motion of the rod and the rigid mass during the impact period. After the cease of impact the motion of the rod is controlled by equation (1), in the same time, since F(t) vanishes, equation (7) declares that the rigid mass moves with a constant velocity.

Finite Element Solutions

The pre-mentioned differential formulation of the equations of motion is equivalent to integral formulation, which requires the application of Lagrange's equation of motion. First, one defines Lagrangean ' L ' by

L = T - [pi] (8)

where 'T' is the kinetic energy and '[pi]' is the potential energy defined by

[pi] = [U.sub.s] - W (9)

[U.sub.s] and W are the strain energy and the work done, respectively, that are given by



W = F(t)u(l, t) (12)

The finite element shape functions ???) (x N link the displacement 'u ' to nodal displacement vector {U} through

u(x, t) = [N]{U} (13)


L = 1/2 [{[??]}.sup.T][M]{[??]}-1/2[{U}.sup.T][K]{U} + [{f(t)}.sup.T]{U} (14)

[M] and [K] are the global mass and stiffness matrices and {f(t)} f is the global force vector. {f(t)} contains only the nodal forces of the last element, [{f(t)}.sub.el], due to the contact force. The rest of the global force vector is full of zeros. [{f(t)}.sub.el] is given by

{f(t)}.sub.el] = F(t)[[N(x = l)].sup.T] (15)

The Lagrange's equation of motion is given by

d/dt ([partial derivative]L/[partial derivative]{[??]}) - [partial derivative]L/[partial derivative]{[??]} = 0 (16)

this leads to

[M]{[??]} + [K]{U} = {f(t)} (17)

Equations (6) and (13) give

F(t) = EA[[partial derivative]N/[partial derivative]x](x = l)][{U}.sub.el] (18)

where [{U}.sub.el] is the nodal displacement vector for the last element, which is the element in contact with the rigid mass.

Equations (17) and (18) are the finite element equations of motion during impact. These equations are applied for both cases of free- free rod and free-fixed rod, see Fig.2. In the case of free-fixed bar both [M] and [K] are positive definite matrices. For the free-free bar, though [M] is positive definite matrix, [K] is positive semi-definite matrix due to the existence of rigid body modes.


Newmark implicit time stepping method (Bathe [19]) is used to express the current velocity [{[??]}.sub.N+1] and acceleration [{[??]}.sub.N+1] in terms of the current displacement [{U}.sub.N+1] and previously determined values of displacement [{U}.sub.N], velocity {[??]}.sub.N], and acceleration [{[??]}.sub.N]. Combining these equations with the equations of motion (17) and (18) yields a system of algebraic equations in terms of [{U}.sub.N + 1] and F[(t).sub.N+1]. Newton-Raphson iterative method (Bathe [19]) is used to solve the resulting equations to find the current displacement and contact force. The displacement and other variables' distributions in the rod at the end of impact serve as the initial conditions for the subsequent free vibrations of the bar, which are governed by the solution of equation (17) while equation (18) is no longer relevant.

Numerical Simulation

Numerical simulations, for a rigid mass collides with a free-free elastic rod and with a free-fixed elastic rod, are presented in this section, see Fig. 2. The rod in both cases is an aluminum rod with a 3 mm x 25 mm cross section, 200 mm length, 70 GN/[m.sup.2] Young's modulus, and 2710 kg/[m.sup.3] mass density. The rigid mass has the same mass as the rod. The mass is moving towards the rod with a velocity of 1 m/s. Fifty elements are used to model the rod in the finite element model. The elements are twonodes and one-dimensional linear elastic elements. The velocity of the created wave is c = 5082.35 m/s and the time for the wave to travel cross the rod from one end to the other end is [tau] = l/c = 3.935 x [10.sup.-5] s.

According to St. Venant's principle, as contact starts the velocity of the contact end becomes immediately equals to the rigid mass velocity and right away a compression wave is created at the contact end and travels across the rod with velocity 'c'. The initial compression stress at the contact end is [[sigma].sub.0] = [v.sub. 0] [square root of E/[rho]] and the stress at that end starts to decrease with time until the reflected wave reaches that end.

For the bar with the other end free, the stress at the free end is always zero, therefore, the traveling compression stress wave is reflected at the free end as a tension wave and whenever that tension wave reaches the contact end at time o, it nullified the stress at that end and contact is terminated. Following the cease of contact, the wave is reflected from the contact end as compression wave and periodic cycles start with a period equals to [tau].

The finite element solutions successfully predict this phenomenon as can be seen in Figs. 3- 7. Figures 4-7 show the dimensionless displacement c/[v.sub.0]l u, velocity v/[v.sub.0], stress -c/[v.sub.0]E [sigma], and contact force 1/[v.sub.0] A [square root of E [rho]] F, respectively, with respect to dimensionless time c/l t. Slight numerical damping is introduced to reduce the oscillations in solutions. Dimensionless analytical solutions are given in [4]. Very good agreement is found between the solutions of the proposed finite element model and the analytical solutions. Figure 3 shows the distribution of the dimensionless stress over the dimensionless length x/l, at equal dimensionless time steps of 0.125. Therefore, the wave propagation can be visualized in that Figure. Figures 3 and 7 show that the arrival of the reflected tension wave into the contact end nullifies the contact force. Therefore, it marks the end of impact. Most of the time, a portion of the rod is in tension while the other portion is in compression as can be seen in Fig. 3. Therefore, the mid-point stress alternate between compression and tension marked by the arrival of the wave at that point, see Fig. 3 and Fig. 6. Though the displacements of the bar ends are continuous, see Fig. 4, the slope of each displacement history suffers discontinuity corresponding to the arrival of the wave at that end, which is reflected in the discontinuity of the velocities, Fig. 5. The time history of velocity in Fig. 5 points up that after the end of impact, the striking mass does not change its original moving direction and the bar starts a continued free vibration. The bar has an average rigid body motion velocity and for each end, the velocity is varying between two limits. The arrival of the wave at each bar end increases the velocity of that end impulsively to its maximum value. The analytical solutions, given by Goldsmith [1], predict that the final dimensionless velocity of the rigid mass to be 0.1353 and the present simulation predicts 0.1469, see Fig. 5. At the arrival of the reflected wave to the contact end the analytical dimensionless stress is 0.1353 and in Fig. 6 the finite element calculates 0.1325.







Using the present finite element simulation, a visualization of the wave motion is illustrated in Fig. 8 for the bar with one fixed end. The figure illustrates the distribution of the dimensionless stress over the dimensionless length x/l, at equal dimensionless time steps of 0.125. It shows that the traveling compression stress wave is reflected at the fixed end as a compression wave, as expected. Since the contact is not terminated yet, the contact end operates as a fixed end and the compression wave is reflected from that end as a compression wave again. Once more, the wave is reflected as a compression wave at the fixed end, but shortly after that the contact is terminated. For that reason, during the contact period the whole rod is under compression all the time, see Fig. 8. Using wave propagation theory, analytical solutions given in [1, 2] predict the same phenomena but without presenting a similar figure. After the end of impact, the subsequent free vibration of the bar has a periodic cycles with a period equals to 2[tau] as can be seen in Figs. 8-11. The free end reflects the wave, as anticipated, as opposite polarity wave. Therefore, during these periods, most of the time a part of the rod is in compression while the other part is in tension, see Fig. 8. The time histories of the bar displacement, velocity, stress, and contact force in dimensionless forms are shown in Figs. 9-12. The figures illustrate that the displacement is continuous while velocity, stress, and contact force suffer discontinuities. At any location in the bar, the discontinuities occur at intervals correspond to the arrival of the waves to that location; see Figs. 8-12. Figs. 8, 11 and 12 confirm that the arrival of the reflected compression wave to the contact end raises the stress at the contact end, and accordingly the contact force, to its maximum value. Next, contact force starts to decrease and impact is terminated when the stress at the contact end vanishes. The analytical solutions given in [1,2] predict the dimensionless duration time, displacement of contact end at separation and its maximum value after separation, and maximum contact force. Both of the analytical results and the corresponding results of the current finite element simulation are tabulated in table1. It has to be noticed that slight numerical damping is introduced to reduce the oscillations in the numerical solutions.


A finite element simulation for the impact of a rigid mass on an elastic rod has been presented in this paper. The impact model utilized St.Venant classical impact model and the two cases of free-free and free-fixed elastic rod have been investigated. As contact established, a wave is initiated at the contact end and starts to propagate through the rod. The wave propagation and the contact force differential equations have been obtain and the finite element discretization of the equations of motion has been developed. Numerical solution procedure has been proposed along the lines of Newmark implicit integration method and Newton-Raphson iterative technique. The current simulation has a valuable advantage over other similar simulations in calculating the contact force according to St. Venant model and also it is not time consuming.

Results show the variation of contact force, displacements, velocities, and stresses with respect to time for both cases of free and fixed far end of the bar. Very good agreement has been found between numerical results and the well-known analytical results. A simulated visualization of the propagation of the stress wave through the bar has been developed. This visualization enhances the understanding of the physical phenomena of impact and wave propagation including the reflection of the wave at free and fixed ends as well as at the contact end. The results demonstrate that the proposed finite element simulation is accurate enough to depend upon for further investigation in wave analysis and simulation.

A cross- sectional area ([m.sup.2])
c wave propagation velocity (m/s)
E Young's modulus (N/[m.sup.2])
F Contact force (N)
{f} global force vector (N)
{f}.sub.el] force vector for element in contact
 with the rigid mass (N)
[K] fglobal stiffness matrix (N/m)
l length (m)
L Lagrangean (J)
m mass of the rod (kg)
[m.sub.0] mass of the rigid mass (kg)
[M] global mass matrix (kg)
[N] finite element shape functions (m/m)
q displacement of the rigid mass (m)
[??] velocity of the rigid mass (m/s)
t time (s)
[t.sub.c] contact period (s)
T kinetic energy (J)
u displacement of the rod at position x (m)
{U} nodal displacement vector (m)
{[??]} velocity vector (m/s)
{[??]} acceleration vector (m/[s.sup.2])
[{U}.sub.el] nodal displacement vector element in
 contact with the rigid mass (m)
[U.sub.s] strain energy (J)
[v.sub.0] initial velocity of the rigid mass (m/s)
W work done (J)
x position in the rod (m)
[epsilon] strain in the rod (m/m)
[pi] potential energy (J)
[rho] density (kg/m3)
[sigma] stress in the rod (N/m2)
[tau] time for the wave to travel cross the
 rod from one end to the other end (s)
[().sub.N] value at the previous time step (N = 0,1,2, ...)
[().sub.N+1] value at the current time step (N = 0,1,2, ...)


[1] Goldsmith, W., 1960, Impact: The Theory and Physical Behaviour of Colliding Solids, Edward Arnold, Ltd., London, UK.

[2] Timoshenko, S.P. and Goodier, J.N., 1970, Theory Of Elasticity, McGraw Hill, New York.

[3] Shi, P., 1997, "Simulation of impact involving an elastic rod. Computer Methods in Applied Mechanics and Engineering," 151, pp. 497-499.

[4] Hu, B., and Eberhard, P., 2001, "Symbolic computation of longitudinal impact waves," Computer Methods in Applied Mechanics and Engineering, 190, pp. 4805-4815.

[5] Hu, B., and Eberhard, P., 2004, "Simulation of Longitudinal Impact Waves Using Time Delayed Systems," Journal of Dynamic Systems, Measurement and Control (Special Issue), 126(3), pp. 644-649.

[6] Keskinen, E., Kuokkala, V-T, Vuoristo, T. and Martikainen, M., 2007, "Multibody wave analysis of axially elastic rod systems," Proc. Instn. Mech. Engrs, Part K: J. Multi-body Dynamics, 221, pp. 417-428.

[7] Maekawa, I., Tanabe, Y. and Suzuki, M., 1988, "Impact stress in a finite rod," JSME International Journal Series I, 31, pp. 554-560.

[8] Hu, B., Schiehlen, W. and Eberhard, P., 2003, "Comparison of analytical and Experimental results for longitudinal impacts elastic rods," Journal of Vibration and Control, 9, pp. 157-174.

[9] Al-Mousawi, M.M., 1986, "On Experimental Studies of longitudinal and flexural wave propagations: An annotated bibliography," Applied Mechanics Reviews, 1986, 39, pp.853-864.

[10] Ueda, K. and Umeda, A., 1993, "Characterization of shock accelerometers using Davies bar and strain-gages," Experimental Mechanics, 33(3), pp. 228-233.

[11] Rusovici, R., 1999, "Modeling of shock wave propagation and attenuation in viscoelastic structure," Ph.D. Dissertation, Virginia polytechnic institute and state university.

[12] Ramirez, H. and Rubio-Gonzalez, C., 2006, "Finite-element simulation of wave propagation and dispersion in Hopkinson bar test," Materials and Design, 27, pp. 36-44.

[13] Allen, D. J., Rule, W. K., and Jones, S. E., 1997, "Optimizing material strength constants numerically extracted from Taylor impact data," Experimental Mechanics, 37(3), pp. 333-338.

[14] El-Saeid Essa, Y., Lopez-Puente, J., and Perez-Castellanos, J. L., 2007, "Numerical simulation and experimental study of a mechanism for Hopkinson bar test interruption," The Journal of Strain Analysis for Engineering Design, 42(3), pp. 163-172.

[15] Seifried, R., Schiehlen, W. and Eberhard, P., 2005, "Numerical and experimental evaluation of the coefficient of restitution for repeated impacts," International Journal of Impact Engineering, 32(1-4), pp. 508-534.

[16] Seifried, R. and Eberhard, P., 2005, "Comparison of numerical and experimental results for impacts," ENOC-2005, Eindhoven, Netherlands, 7-12 August, pp. 399-408.

[17] Krawczuk, M., 2002, "Application of spectral beam finite element with crack and iterative search technique for damage detection," Finite Elements in Analysis and Design, 38(6), pp. 537-548.

[18] Krawczuk, M., Palacz, M and Ostachowicz, W., 2004, "Wave propagation in plate structures for crack detection," Finite Elements in Analysis and Design, 40(9-10), pp. 991-1004.

[19] Bathe, K-J., 1996, Finite Element Procedures, Prentice Hall, Upper Saddle River.

Biography for Dr. Hesham A. Elkaranshawy

He received his BSc. and MSc. degrees in Mechanical Engineering and in Engineering Mathematics from Alexandria University, Alexandria, Egypt, in 1983 and 1988 respectively, and his PhD degree in Mechanical Engineering from McMaster University, Hamilton, Ontario, Canada, in 1995. He worked in CAE Electronics, Montreal, Canada, during 1995-1996 as a finite element specialist for its part in the world space station project. He joined Department of Engineering Mathematics, Engineering Mechanics division, Alexandria University as an assistant professor in 1996. In 1998 he established Dynamica Engineering Consultants, Giza, Egypt. Currently, he is an assistant professor at Department of Mechanical Engineering, Yanbu Industrial College, Royal Commission for Jubail and Yanbu, Saudi Arabia. His research interests include multibody dynamics, dynamic modeling of robotic systems, contact-impact dynamics, finite element modeling, mechanical vibrations, wave propagation and fault detection.

Hesham A. Elkaranshawy (1) and Nasser S. Bajaba (2)

(1) Dept. of Mechanical Engg. Technology, Yanbu Industrial College, Yanbu, Saudi Arabia, P.O.: 30436 e-mail:

(2) Dept. of Mechanical Technology, Jeddah College of Technology, Jeddah, Saudi Arabia e-mail:
Table 1: Comparison between the analytical and the finite element
results (free-fixed rod)

 Analytical Proposed Finite
Dimensions Values Results Element Results

Duration Time 3.068 3.075
Displacement of Contact
 end at separation 0.375 0.371
Maximum Value of
 displacement of contact
 end after 0.576 0.576
Maximum Contact Force 2.135 2.045
COPYRIGHT 2008 Research India Publications
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2008 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Elkaranshawy, Hesham A.; Bajaba, Nasser S.
Publication:International Journal of Applied Engineering Research
Article Type:Report
Geographic Code:1USA
Date:Jun 1, 2008
Previous Article:Estimation of solar irradiance for PV--ECS based distributed power generation.
Next Article:Analytic and FEM evaluation of power density for various types of double-sided axial flux slotted PM motors.

Related Articles
Using finite element analysis in designing roadside hardware.
Elastomer Biaxial Characterization Using Bubble Inflation Technique. II: Numerical Investigation of Some Constitutive Models.
Finite element modeling of polymer hot embossing using a glass-rubber finite strain constitutive model.
A review on cold forging die design and die design process.
Micromechanical analysis of hybrid discontinuous fiber reinforced composite for longitudinal loading.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters