A study of axial impact of composite rods using SPH approach

A computational simulation of axial impact between composite rods is presented in this paper. The smoothed particle hydrodynamic (SPH) method is coded and implemented for this study, in which the property and strength for different materials are included. The SPH formulation for impact between rods as well as its computatinal procedure is also expounded. The present approach is used to study impact between composites rods parametrically. The effects of material properties and configurations of composite rod on impact response are then examined in detail. In addition, the present SPH results are compared with those from theoretical analysis.


Introduction
Smooth Particle Hydrodynamics (SPH) is a technique for problem solving in Computational Continuum Dynamics [1].Compared with the finite element method (FEM), it's meshless and thus can effectively avoid mesh crash problem and remesh process, which is especially suitable for the problem involving large deformations.This method also adopts a Largrangian frame of reference, which is direct and simple to manipulate.Because of its simplicity and versatile, the SPH method is receiving much attention recently.
Smooth Particle Hydrodynamics (SPH) was initially developed in astrophysics by Lucy [2], Monaghan and Gingold in 1977 [3].Soon the SPH method was applied to liquid problem and liquid-like problems such as hypervelocity impact.It was then extended to treat the dynamic response of solids by Libersky et al. in 1991 [4], in which Mie-Gruneisen equation as the equation of state for solid and formula for deviatoric stress component was adopted.Libersky [4] verified his im-plementation with Noh problem.Johnson [5] used cylinder impact problem to study artificial viscosity effect.Libersky [6,7] compared cylinder 2D coordinate system with 3D implementation using cylinder impact.Dyka et al. [8] presented stress point method through a simple elastic bar.Swegle proposed an impact problem involving dissimilar materials as a test for SPH.In Libersky's paper [1], he used Swegle's 1D problem to verify his boundary treatment.
Although the SPH method has been widely applied in fluid and solid dynamic analysis, the simulation of impact between rods using SPH approach has not been fruitful, especially for rod made of composite materials.This paper address this and present an feasible SPH approach to simulate axial impact between composite rods.The merits of SPH method is exploited for this study and also compared with those from theoretical analysis.The present SPH approach is used to examine effects of material properties and configurations of composite rod on impact response; the compuational results are illustrated and discussed in detail.

Problem description
Consider three cases of axial rod impact as shown in Table 1.Both the striking rod and the target rod has a diameter of 2 mm.In Case 1, the striking rod, made of copper, has a length of 10 mm.It impacts the target rod at an initial velocity of 20 m/s.The target, rod with a length of 40 mm, is static and free along the x-direction initially.Four types of the materials for the target rod are considered in the first case.They are copper, steel, aluminum and polymethyl methacrylate (PMMA) respectively.
In Case 2, the geometry, material and initial condition of striking rod is the same as the striking rod in Case 1. Unlike the axial impact between the rods in Case 1, the striking rod will impact the rigid wall.In Case 3, a pad of PMMA is attached to the head of the striking copper rod with a length of 10 mm.The target copper rod, with a length of 40 mm, is also static and free along the x-direction initially.The length of the PMMA pad may vary whereas the length of the striking rod is kept unchanged in this case.
Table 2 shows the properties of materials used in the above three case studies.ρ 0 is the initial density, C 0 is the elastic wave speed, S is the experimental data from Hugoniot curve, Γ is the Gruneisen constant and Y 0 is the yield stress.PMMA is used for its merit of low shock impedance (ρ 0 C 0 ) and such impedance is very close to the sensor used to detect the stress.

Governing equation
For one dimensional continuous media mechanical problem, the conservative equation for mass, momentum and internal energy are as follows.
Where x is the Lagrangean coordinate, t is time, v is the particle velocity and ρ is the density; σ and e denote the normal stress in the x-direction and unit internal energy, respectively.
In order to solve the problem, the equation of state is needed.
In which p is pressure, p H (η) can be expressed as Where η is relative compression, defined as: The relation between stress and pressure can be expressed as: For one-dimensional question, there is σ = −p.The stress can be determined from Gruneisen form of equation of state, which is written as: Where the H subscript refers to values on the Hugoniot curve, v is volume and Γ(v) is the Gruneisen parameter defined as which can be interpreted as a statement that the pressure and energy at any particular volume can be related to these quantities on the Hugoniot through some average value of the Gruneisen parameter.
For most solids and many liquids over a wide range of pressure (within one phase) experiments have shown the existence of an empirical linear shock velocityparticle velocity form of equation of state (c = C 0 + Sv).Solving this equation, in conjunction with the conservation equations (mass, momentum and energy) for p and e, we can get: Substitute these into Eq.( 9), we can get the Mie-Gruneisen equation of state Eq.(4).
Equations ( 1)-( 3), ( 8) in conjunction with Eq. ( 9) reveal the relationahip between stress, density, velocity, energy as well as material perpoties of the rods, which provides a theoretical base for the study of the axial impact between composite rods.

SPH discretization of governing equations
SPH method is used for the discretization of governing equations.The basic formula for SPH discretization can be expressed as: This is the most important formula in SPH method, where W i (x − x i , h i ) is the kernel function.The key concept is that W has a compact support, which means only a limited number of neighbor points are involved in the evaluating of one point.The widely used kernel function W is the spline function.For 1D problem, it is as follows: for q 2 where q = s/h.In order to discretize the governing equations, we need the first order derivative of the functions which also can be approximated by the derivative of SPH interpolating of the function.It means that only derivative of kernel function is needed, that is: For simplicity, following notations are used throughout the deduction.
From the basic form, let f = ρ, it's natural to get the following representation of density The density at any point in space is then obtained by summing up the contributions from all particles at that point.SPH derives its name from this interpretation.Equation (18) requires only particle coordinates and masses to compute the density and automatically satisfies the continuity equation, provided that the particle masses are constant.
Actually, we can also discretize conservative Eq. (1) of mass directly to get another discretization of mass conservation.It's reported by Libersky [6] that the density-by-summation form can provide greater accuracy and robustness.So here we will not discretize the mass conservation equation.
For the momentum conservative equation, using the SPH discretization Eq. ( 13) and integration by part, we can get the discretized conservative form of momentum equation: With the same approach, we can get internal energy conservation: The above forms of discretization are widely adopted, mainly due to they are symmetric about i and j, which help to keep the conservation.

Artificial viscosity
The artificial viscous stress ij derived by Monaghan and Gingold is used in this study [9], which is as follows: Where The average of elastic wave speed and average of density Eq. ( 24) are used to ensure the conservative condition.Artificial viscosity changes abruptly around the discontinuity.

Computational procedure
Both striking and target rods are discretized into a row of smoothed particles evenly distributed.Initially particles representing the striking rod have a velocity of 20 m/s.By assuming the cross-section of the rod to be unit, we can get the mass of each particle.

Boundary and interface
For the problem of axial impact between composite rods, three types of boundary conditions are considered.They are rigid wall, free boundary and stress-free boundary.The rigid wall boundary is handled by the ghost particle method.The ghost particle method actually converts the problem that the striking rod impacts the rigid wall into a symmetric impact.
For free boundary stress-free boundary, density deficiency is the main problem.By using the kernel estimation of unity centered at the boundary particle, Liber-sky [1] derived formula for the density of a boundary particle i with the kernel summation corrected for boundary deficiency The above algorthm is used in this study, and therfore the density at the boundary will be evaluated only with the mass density and position of the internal particles, the so-called density defficiency will be overcomed.
The interface between different materials was dealt with by a smoothing method to evaluate the density.The procedure of the method is that an initial smooth on the initial density is done before the main computing cycle and the smoothed density as the initial density in the main computing cycle is taken.Therefore, the density will not change due to the abrupt change of density around the interface and thus the pressure will not get the unphysical change from the interface.

Time stepping
Time integration is another important aspect in the computational procedure.Currently, we use leapfrog method to do time integration.Leapfrog method was a conventional two-step method.This method can ensure second order accuracy.The advantage of the method is its low memory storage and efficiency (one force evaluation per step).And the method is dependent of the Courant-Friedrichs-Lewy condition for stability, which typically results in that time-steps are proportional to smoothing lengths.When particles' clustering occurs, the smoothing lengths become very small, and then this method will become prohibitive.But to our solid response problem, this will never happen.
Following is an illustration of the two time-steps.First step: Second step:

Correlative studies
First, a problem of a short rod impacting a long rod is simulated.Both short and long rods are composed of copper.This case study is conducted to verify the implemented SPH code and also to check the suitability of the SPH approach for the problem of the axial impact between rods.The results obtained from the SPH approach were compared with those from theoretical analysis.
For the rod impact, the relation between wave speed, particle velocity and stress can be found from the water hammer equation [10]: Where [σ] is the change of the stress, c is the wave velocity which takes account of the direction of propagation, i.e. c = +c if the wave propagates in the positive x-direction and c = −c if it goes in the negative x-direction.[v] is the change in particle velocity, taking into consideration the sign of the velocity change.
At the contact surface during impact, the contact force between the rods is common, so there is the common stress σ c since the cross-section areas of the striking rod and target rod are the same.From the contact condition, there should also be common particle velocity v c in the portions of the rods where stress waves initiated.For the case where both striking rod and target rod are made of copper, substitute the condition for both striking rod and target rod into Eq.(30) we can get, Where v s is the initial velocity of the striking rod, and v t is the initial velocity of the target rod, here we assume v t = 0. Solving Eqs (31)-(32), we can get, v c = v s /2 and σ c = −ρcv s /2.Where σ c is the common stress spread from the contact surface, v c is the common particle velocity.
Assume the contact surface at a position of x = x 0 and the target rod has a length of l.The stress wave propagates in the target rod, when −x 0 < x − ct < x 0 and 0 < t < l/c with σ = σ c and v = v c .Condition for t < l/c means this shock pulse response is applicable only before the shock front hits another free end of the target rod.
Figure 1 illustrates the stress distributions at time t = 7.7 µs from both theoretical analysis and the SPH approach.It is seen that the length of the pulse is about twice as that of the striking rod. Figure 2 shows the particle velocity distributions at a specific time t = 12.8 µs from both theoretical analysis and the SPH approach.At this instant, the stress is zero along both striking rod and the target rod.The particle velocity around the rear part of the target rod goes up to the initial velocity of striking rod and the length also equal to that of that striking rod.This result coincides with the superposition method, which was widely used.The slope around the abrupt change can be explained by the fact that SPH kernel approximation can not simulate the abrupt change.
From Figs 1 and 2, it is seen that the results from the SPH approach match reasonably up those from the theoretical analysis.This indicates that the SPH approach yields acceptable results for the problem of axial rod impacts.

Effects of material properties on resulting stress
The computational simulation of axial rod impact is conducted using the SPH approach.The short copper rod impacts the long target rod and rigid wall at a velocity of 20 m/s respectively.The materials of long target rod are assumed to be copper, steel, aluminum and PMMA respectively as illustrates in Table 1. Figure 3 shows the resulting stress history at the center point of the striking rod (x = 5 mm) for the five different targets: rigid wall, long rod made of copper, steel, aluminum and PMMA respectively.From Fig. 3 it is observed that the stress in the striking rod varies with the material properties of the target.The maximum stress value for impacting on the copper rod is reduced to 50% of that for impacting on the rigid wall.And further more, the maximum value for impacting on PMMA rod is only 8% of that for impacting on the rigid wall.
Figure 4 shows the resulting stress history at x = 15 mm of the target copper, steel, aluminum and PMMA rods respectively.For the copper target, the stress is reduced most rapidly to zero.For the steel and aluminum target, the stress is also reduced quickly with a very low magnitude of stress remnant.For the PMMA target, the stress response is almost kept unchanged over the most time duration.
Figure 5 depicts a particle velocity distribution at time t = 5.1 µs, when the first stress pulse moves out the striking rod.The difference between the shock front is due to the different wave speed of the target  materials.It can be seen from Fig. 5 that after this moment, the copper striking rod gains different speeds corresponding to each target rod with different materials.When the target rod is made of same material, striking rod totally stopped, which means all the energy and the momentum goes into the target rod.On the other hand, when the copper rod impacts a PMMA rod, it retains a speed about 85% of the initial striking speed; the remaining momentum will be transferred to the target gradually.And the impact duration for the impact between the copper and PMMA rods will be much longer than the duration for the impact between the same copper rods.
From Figs 3-5, it can be found that the shock impedance of the target has a significant influence on maximum stress in the striking rod.The stress in the striking rod may decrease with the decrease of the shock impedance of target.It is also observed that the worst condition is the striking rod impacting a rigid wall, in which the shock impedance tends to infinity.

Effects of the PMMA pad on attenuation of impact response
The present SPH approach was also used to study the effects of the PMMA pad on resulting stress.As illustrated in case 3 (See Table 1), the PMMA pad is attached to the head of copper striker.The length of the PMMA pad rod L p is 1/2, 1, and 2 times length of the striking rod respectively.
Figure 6 shows the resulting stress histories of the copper striker at x = 5 mm for three lengths of the PMMA pad rod. Figure 7 is the resulting stress histories of the PMMA pad with three different lengths.A comparison between the stress curves is made for both cases of the copper striker attached the PMMA pad rod (L p = 2L s ) shown in Fig. 6 and the copper striker without any attachments shown in Fig. 1.It is observed that the maximum stress is reduced by 74% when the copper striker is attached a PMMA rod.This indicates that the PMMA pad has a significant effect on the attenuation of impact response.
From Figs 6 and 7, it is also observed that the length of the PMMA pad influences the magnitude of the stress response.When the length of the PMMA is less than the length of the copper striker, i.e.L p L s , the stress response decrease rapidly as the length of the PMMA pad increase.When the length of the PMMA rod is larger than the length of the copper striker, i.e.L p > L s , the effect of the length of the PMMA pad is insignificant.
These results show that a pad with a low shock impedance material can be used to attenuate the maximum stress in the striker induced by the impact.

Discussion
As the compressive wave reflects and becomes tensile wave, the free end will keep stress-free, so the method can simulate this process exactly.Unlike theoretical analysis, the SPH method is not required to specify that boundary should be stress-free, the SPH approach makes the stress zero automatically.
After the first computing cycle, the density would be smoothed by the mass conservation formula, and according to equation of state, a pressure and particle velocity would be induced since density is changed,which is an unphysical procedure.This is the main problem of the interface.Since we changed the initial density, which is also used in the equation of state, according to the first round summation computing for density, then there is no unphysical pressure and particle velocity produced.The following density by summation proce-  dure will not change the already smoothed density distribution.It is noticed that when a shock pulse pass the smoothed interface, the unphysical pressure still exists, but the magnitude is much smaller.The smear of the shock and release front is caused by both the kernel estimate and the equation of state.
A longer smoothing length, e.g.including more neighbour particles, can effectively control tensile instability, but this will alternatively smear the abrupt change of shock front.It's the defect of SPH interpolation itself, because the value is computed over the neighbor particles, more particle involved, more vulnerable to disrupt the steep shock front.

Conclusion
The SPH approach for the simulation of axial impact between composite rods was presented in this paper.The smoothed particle hydrodynamic (SPH) method was coded and implemented with the property and strength of different materials for the study of axial impact between composite rods.The computational results show that the material properties of the striking rod have significant effects on the resulting stress response.The PMMA pad can be used to attenuate impact response and further more the length of the pad can influence the attenuation effects.The present SPH results were also compared with those from theoretical analysis, which indicates that the SPH approach can be used for the rod impact problem.

Table 1
Geometry of striking and target rods and coordinate system used for the three cases