Numerical simulation of axially impacted elastic bar

This numerical analysis considers a model of an elastic cylindrical bar impacted by rigid mass, where the accuracy of the model is benchmarked by a closed-form solution developed for this impact process. The analytical solution provides theoretical background and understanding, however, practical problem is not always restricted to a simple one-dimensional geometry and boundary conditions of a purely elastic material behaviour. Hence, the current work numerically simulates the impacted elastic bar, then explores parameters which are not considered in its theoretical model. A two-dimensional (2D) solid axisymmetric system is considered, in which rigid impactor is assigned with initial velocity corresponding to certain drop height and a 2D contact type is defined between interacting bodies. Both elastic and time-dependent viscoelastic material models are considered in this study. The simulation results reveal time intervals gradually increase in every sequential intervals, while relatively small discrepancies are recorded for peak load and pulse width outputs. The parameters of impactor’s mass, drop height and structural stiffness are varied; showing how these parameters individually affect the resulting force response at end struck. The models then evaluate time-dependent viscoelastic material model; showing stiffer response of the resulting force in comparison to models assigned with long-term elastic moduli. Derived elastic modulus-output variable relations show comparable mathematical forms to corresponding equations formulated analytically.


INTRODUCTION
This numerical work simulates an impact problem based on Timoshenko and Goodier [1], which involves a restrained elastic bar subjected to rigid mass as shown in figure 1. The simulation is undertaken using a finite element (FE) commercial package, LS-DYNA (smp d R7.1.1, LSTC, Livermore, California). The available analytical model is considered as the bar is treated as purely elastic, then it is extended to more complex time-dependent elastic behaviour. The study also evaluates relations between several parameters and corresponding output variables.
The dynamic properties of trabecular bone were characterised by Shim et al. [2]: the time-dependent viscoelastic properties using one Maxwell element in parallel with a nonlinear Newtonian dashpot constructing a nonlinear viscoelasticity model. The one-dimensional (1D) model was extended to a threedimensional (3D) model by incorporating 3D strain tensor into 1D strain variable in the non-linear strain rate dependency constitutive model.
Time-dependent properties of bovine trabecular bone samples were characterised in Manda et al. [3] using uniaxial creep tests. The research utilised the creep compliance function, which is easily convertible to the relaxation modulus function, E(t) by using numerical interconversion methods [4]. IOP Publishing doi: 10.1088/1757-899X/1244/1/012015 2 As a result, three-term Prony series of compliance-relaxation functions were derived describing linear viscoelastic response of the bone specimens subjected to compressive creep load, which are related to their bone volume fraction (BV/TV).

PROBLEM DEFINITION
Consider a stationary elastic cylinder with diameter, 2r, and length, l, subjected to impact due to an incoming rigid impactor at one end (called the end struck), while the other end is restrained as shown in figure 1. Numerical solutions are important as not all problems are elastic, have a simple one-dimensional geometry and boundary conditions as in theoretical model by Timoshenko and Goodier [1] and Rosli [5]; consequently, closed-form solutions are not possible. The closed-form solutions, however, help in evaluating and benchmarking the accuracy of numerical methods, such as this attempted work. Explicit code of LS-DYNA is utilised for this work, which uses explicit time integration based on central difference method (CDM) to obtain nodal displacements from their velocities and accelerations in every time step. The time step is calculated by the programme at every stage in simulation by dividing the smallest characteristic length with the continuum wave speed [6]; indicating the time required for wave speed to traverse the smallest element. The current study is limited to elastic and viscoelastic material models without modelling their post-elastic behaviour.

NUMERICAL MODEL
FE analysis requires four basic main inputs: the geometry with corresponding FE mesh; boundary conditions defined for the models; load application and the material behaviour of the models. This section discusses the definitions for each of these required inputs in the current study.

Discretisation of axisymmetric system
Impact mechanics involves more than one body in the collision [7], which implies the necessity of defining a minimum of two geometrical models. Planar symmetry of the cylindrical shape in the defined problem as shown in Fig. 2 is exploited. Two-dimensional (2D) solid models in LS-DYNA are based on integral difference scheme, which defines the geometries in global XY plane. In axisymmetric case, radial direction corresponds to X-axis, while the axis of symmetry lies on the Y-axis. LS-DYNA offers two options of axisymmetric elements, namely Petrov-Galerkin and Galerkin finite element approaches, generally referred  The volume-weighted Galerkin method is chosen in this work since there is no possibility of a hydrodynamic scenario arising in the problems being considered. The bar is discretised into 4-noded axisymmetric 10 × 40 elements to provide almost square 2D elements, with similar 1:2 diameter-to-length ratio (figure 2b). The rest of this primitive geometry are theoretically revolved 2π radians around the axis of symmetry to create a virtual cylindrical model of the bar as shown in figure 2.
Accordingly, the impactor is also assumed to have axisymmetric shape with arbitrary radius-length dimension, while ensuring the impactor's mass is 2.5 kg simulating the experimental drop test or around 26.6 g for a light-weight dropped hammer [5,6]. Figure 2(a) shows the impactor as coarsely meshed in order to reduce the computational cost, while refinement is applied at the bottom elements which involve with contacting the end struck of the bar. This mesh refinement is required to ensure each slave node (on the impactor) are properly projected to its master segment (the bar struck). Node-collapsing method is applied to some three-noded elements in this region of refined mesh.

Boundary conditions
By the definition of axisymmetry in LS-DYNA, the global Y-axis represents the axis of symmetry, which for the problem considered automatically constrains the nodes on this axis in the global X-translation (radial). In addition, all bottom nodes are restrained in vertical direction as shown in figure 2(b) which is sufficient and consistent with the closed-form model in figure 1.

Loading conditions: Contact and initial velocity
Contact definition is required to establish interaction between the contacting surfaces of impactor and bar. The contact definition allows the incoming slave nodes of the impactor to recognise the top of the bar (which contains master segments). The mesh refinement size between the two contacting models is comparable. Standard penalty-based contact is employed for this purpose, which uses numerical springs to prevent penetration of each slave nodes through master segments as well as to transfer loads between contacting parts [6].
The 2D Automatic Surface to Surface contact type is employed which uses penalty forces to prevent penetration between external faces of 2D continuum elements. In addition, Force Transducer option is used in conjunction with the previous contact card to measure forces generated by that contact definition.
In the current axisymmetric model, the force output is per unit radian. Normal direction of both slave IOP Publishing doi:10.1088/1757-899X/1244/1/012015 4 and master surfaces are determined automatically by the programme, in which the initial model set up already has orthogonal contact surfaces prior to the impact [6].
The hammer is assigned with corresponding initial velocity from drop height with only small gap is introduced between the models in order to avoid initial penetration leading to instability in the calculation of the contact forces [6]. The gap size is generally balanced between minimising initial computational cost (caused by a large gap) and avoiding initial penetration (due to gap being too small) between contacting solids. The best decision depends on the pre-calculation of the defined time step of the smallest element, the initial distance between the first expected contacting elements and the assigned nodal velocity.

Material models
The bar's material is defined as elastic, which is similar to the closed-form solution in [1,5], or viscoelastic. In both cases, the bar has homogeneous and isotropic mechanical properties. The impactor is defined as rigid material which matches the problem definition and justifies the irregular mesh of the impactor near the impacted region. In addition, this material model is a cost-effective option because its model does not store history variables due to its undeformed body restriction.
Even though the rigid definition itself implies infinite stiffness, values of both Young's modulus, and Poisson's ratio need to be assigned to determine the bulk modulus for calculating the contact stiffness. In addition, the explicit time integration requires the mass density of the impactor in order to construct mass matrix in the CDM solution for the nodal displacements, even though the wave propagation is not calculated for the rigid material.

COMPARISON WITH THE CLOSED-FORM SOLUTION
The force response of a cylindrical FE model (10.6 mm diameter and 12 mm length) impacted by a light-weight 26.6-g hammer is first compared with the closed-form (CF) solution in terms of the resulting peak load and pulse width. The solid axisymmetric bar is converted to a one-dimensional model by assuming the Poisson's ratio to be zero, i.e. there is no transverse response in the 2D solid element. The light hammer is chosen to increase the mass ratio, which in turn would result in only seven stress intervals in the pulse response, allowing easy comparison to be made.  [5]. The effect of a non-zero FE model is also shown. Figure 3 compares the FE response (zero Poisson's ratio) to that from the closed-form solution. It can be seen that there are seven force intervals in both FE and closed-form responses, which they are similar in the form of bell-like shape. The sharp peak-and-trough of vertical shape in every interval due to the instantaneous rise of force by CF solution [1,5] is not recreated by the numerical simulation, nonetheless it is still fairly visible.
It is apparent that the closed-form solution returns consistent time interval, following strictly the governing equation of interval time [1,5], whereas this time from FE solution gradually increases in every IOP Publishing doi:10.1088/1757-899X/1244/1/012015 5 sequential interval. This discrepancy is explainable by understanding the ideal scenario in the closed-form solution in which the supposed deformation of bar is neglected in the analysis. On the other hand, the FE model experiences numerical damping which is a default setting for all contact options in the LS-DYNA [6] causing longer time for the travelling stress wave. Table 1 shows comparison of the two main outputs under investigation, namely peak load ( !"# ) and pulse width ( $%&'( ) between the two solutions. The FE method records 11.43% lower peak load, while its pulse width is 6.72% larger in comparison to its benchmark results by the CF solution. The slightly high discrepancy recorded by this FE method is explicable by reevaluating the element size of elastic bar which has 1:2 ratio of rectangular shape. It will be shown in the subsequent sections when the elements have 1:1 ratio, the outputs improve significantly as squarelike 2D elements give more accurate results. It is worth mentioning that in both groups of bar models, the total number of elements are similar viz. 400 which are considered as sufficiently refined. The effect of Poisson's ratio (ν = 0.33) is shown in figure 3, in which higher frequency distortion in its overall resulting pulse pattern is observed causing slightly unclear pattern of seven intervals. This happens perhaps due to the effect of wave bouncing from the sides. Nonetheless, the following peak load and pulse width values are not significantly affected, such that only 1.5% and 2.2% differences with the one-dimensional FE model. It can be concluded that assigning a reasonable value of Poisson's ratio does not drastically alter the response, in particular the output variables under interest, which are !"# and $%&'( .

EFFECT OF DIFFERENT PARAMETERS ON THE RESULTING ELASTIC PULSE
The effect of varying several parameters: mass of impactor; drop height; and structural stiffness; towards the resulting force at end struck by using numerical simulation is considered. This elastic study of axial impact on trabecular bone samples is based on report by Rosli [5]; some of the parameters are not varied i.e. density: = 1.31 × 10 )* kg/mm + ; bar's length, = 21 mm; and Poisson's ratio, = 0.33. The axisymmetric model is discretised using 10 × 40 elements in radial and axial directions respectively, which is almost 1:1 of bar's diameter-to-length ratio or even slightly less for models with smaller radius. The analysed output is the summation of eleven (11) nodal forces histories of the nodes at the top of axisymmetric bar model, multiplied by 2π to obtain the total force acting on the surface of end struck.

Impactor's mass and drop height
Two values of impactor's mass are chosen, which are 1.807 and 2.5 kg, while the two drop heights are 100 and 200 mm, making up four (4) FE models for this analysis. The models have geometry of 10.6 mm of diameter and 21 mm of length, with elastic modulus is assigned to be 531 MPa as reported by Xie et al. [8] and falls within ranges reported by many studies including [9][10][11][12][13][14] for trabecular bone. The drop height affects only the peak load, while pulse width is independent of this input parameter as shown in figure 4; this finding by numerical method is consistent with theoretical parametric analysis in Rosli [5]; it was also shown in the study the peak load has direct proportionality to the square-root of drop height, h. As seen in figure 4 , this is satisfied by the numerical solution. The pulse width is also almost identical for impactors of the same mass dropped from different heights. Figure 4 also shows that a higher drop height results in a fractionally smaller pulse width. A larger mass of the impactor leads to larger peak loads and also pulse width, which is consistent to the analysis of the closed-form solution in Rosli [5]. In addition, it is observed in figure 4 that with the increase in ,!$"-./0 , the force responses appear to be almost similar initially, then the load for higher mass increases at a substantially higher rate and as expected comprises of larger number of stress intervals.

Structural stiffness
The structural stiffness of a simple column under compressive load is a function of its material (Young's modulus) and dimensions (cross-sectional area and length). The pulse exhibits stiffer response when the combination of higher peak load and lower pulse width is obtained [5].
In this parametric study, 2.5 kg impactor is dropped from 100 mm height. The parametric study with respect to length was examined in the closed-form solution [5], therefore only elastic modulus and diameter (representing area) are varied in this analysis. Similar elastic modulus of 531 MPa is chosen, with addition of arbitrary materially-stiffer models of 1 GPa. The diameters are either 10.6 mm [5] or 8.6 mm, in which a smaller area is expected to exhibit a less stiff behaviour. The resulting force responses at the end struck are shown in figure 5 for four varying-stiffness models. The stiffest response is shown by the model with highest elastic modulus and cross-sectional area, i.e E = 1 GPa, radius, R = 5.3 mm; followed by models with lower structural stiffness contributors. These demonstrate the influence of elastic modulus and cross-sectional diameter on the peak load and the pulse width. The peak load is seen to be directly proportional to the elastic modulus as would be expected [5]. The peak load is also directly proportional to the cross-sectional area. The effect on pulse width is simply the opposite to that for the peak load i.e. pulse width decreases with increase in elastic modulus and radius.

IMPACT RESPONSE OF MODELS WITH TIME-DEPENDENT MATERIAL BEHAVIOUR
The effect of employing a time-dependent material model, viz. viscoelasticity (VE) is analysed in comparison to simple time-independent linear elastic model (abbreviated as TI-E) assigned to the impacted axisymmetric bar. In VE material model, the relaxation modulus can be defined using the Prony series [4] as (1) where Eeqm is the equilibrium modulus and 2 and 2 are the relaxation strengths and relaxation times respectively. A three-term Prony series ( = 3) was shown to adequately represent the relaxation modulus of trabecular bone with different volume fractions, BV/TV [3]. The values determined in this above cited study are used to consider the influence of time-dependent material properties on impact response.

Implementing viscoelasticity in LS-DYNA
In TI-E model, single value of elastic modulus, E is taken as (∞) or (1! from equation (1); a longterm modulus representing the material resistance value when the load stabilises hence treated as in statically-applied manner in equilibrium. On the other hand, the viscoelastic material model in LS-DYNA, e.g. Mat076: General Viscoelastic requires elastic bulk modulus, K. Similar to the equilibrium modulus in TI-E, K in VE is defined as a long-term modulus. The Poisson's ratio was defined as 0.33. In order to exhibit relaxation phase in the VE model, shear relaxation moduli, 2 and decay constants, 2 for each term in Prony series are required. Hence, the general expression of the Prony series in LS-DYNA is (2) Comparing these two relaxation expressions in equations (1-2), the relaxation strengths and times parameters in equation (1) provided by Manda et al. [3] need to be converted onto shear moduli and decay constants in equation (2) respectively. Shear moduli are obtained from E and ν assuming isotropic behaviour while the decay constants (βi) are simply inverse of the relaxation times (ρi) in which = 1, 2, 3 for three terms of Prony series. In the case of three-term Prony series, the General Viscoelastic material model requires four input terms, in which the extra term is (1! with zero decay constant, fully expressed as Six (6) bone volume fraction models, i.e. BV/TV=19, 26, 33, 39, 42 and 43% are selectively chosen from Manda et al. [3] on the basis of signifying direct correlation between material's stiffness and the reported BV/TV. The required input parameters are given in table 2 by using moduli conversions into equation (3). All these conversions are applicable in this study assuming the material as having isotropic behaviour.

Force response
The bar is discretised into 4-noded axisymmetric 10 × 40 solid elements impacted by a 2.5-kg rigid impactor. The length, l of the bar is 21 mm and its diameter, 2r is 10.6 mm. The impactor was assigned with axial velocity of 0.99 m/s, which is equivalent when it is dropped from 50 mm height. In order to terminate the simulation, the pulse duration is estimated [1,5] as then the termination time was set to be slightly longer than the estimated pulse width in equation ( 4 ).
Two main outputs are investigated from the resulting force response at the end struck of axisymmetric bar, namely the peak load as the highest magnitude of compressive load experienced by the end struck and the pulse width, which is the contact duration of the impact incidence. Fig. 6 shows graphical elastic and IOP Publishing doi:10.1088/1757-899X/1244/1/012015 9 viscoelastic responses of three different BV/TVs, clearly depicts that viscoelastic models have stiffer response than their TI-E counterparts. Figure 6. Effect of varying time-independent and viscoelastic properties towards peak load and pulse width.
The key numerical results for all models are tabulated in table 3 with E(t) represents Eeqm and E(0) for TI-E and VE respectively. The modulus function E(t) represents well the contribution of BV/TV to bone sample's stiffness for both material models. In addition, the stiffness of VE models, which is E(0) is consistently higher than TI-E (Eeqm).  Maximum compressive stresses from table 3 recorded higher values than these reported strengths, nonetheless they were proven experimentally to be realistic in drop tests [5].

Modulus-to-output variable relations
The parametric study in [2] points to the requirement of using linear direct best fit regression for relating √ to the peak load i.e. !"# ∝ √ . On the other hand, its effect on the pulse width is inversely proportional or negative power law i.e. $%&'( ∝ 1 √ ⁄ . The mathematical relations obtained from this numerical analysis is shown in figure 7. The output variable-elastic modulus relations in figure 7 show very close resemblance to the mathematical relations derived from the closed-form solution provided by Rosli [5] i.e. similar regression models return almost identical constants. Hence, it can be concluded that the current FE models well replicate the theoretical solutions, hence they are benchmarked and could be used for more complex axial bar impact problems.

SUMMARY
This numerical work replicates an axial bar subjected to rigid mass impactor problem as described in previous studies. Sufficiently-refined volume-weighted Galerkin axisymmetric elements model is used which is available in LS-DYNA as a cost-effective geometrical model, while linear elastic and timedependent viscoelastic material models were evaluated. The study is successfully validated by theoretical solutions in which comparable results were obtained when analysing: 1) the effect of Poisson's ratio (difference with one-dimensional model is less than 3%); 2) parametric study on impactor's mass, drop height and structural stiffness; and 3) output variable relations (peak load and pulse width) with square-root of elastic modulus. The output variable-input parameter relations are useful to explain empirical mathematical relationships. Hence, the FE models developed in this study are successfully benchmarked and can be further developed for incorporating more complex inputs from axial bar impact problems.

ACKNOWLEDGEMENT
The first author is indebted to the Ministry of Education, Malaysia for the award of PhD sponsorship.