Abstract
There are a few numerical simulation methods available for impact problems. However, most numerical results are not validated experimentally. The goal of this paper is to examine how well the simulation results correspond to the physical reality. In this work, normal and oblique impacts of a hemispherical-tip rod on a square plate are investigated both numerically and experimentally. In the numerical approach, finite element method is used to discretize the contact bodies to describe the deformation precisely combined with the floating reference frame method to describe the rigid motion. In the experimental study, strain gauges and Laser Doppler Vibrometers are employed to measure the high-frequency impact responses. Detailed comparative studies between numerical and experimental results are performed. In the case of normal impact, great attention is given to investigate the influence of finite element mesh size on the simulation accuracy and a “Prediction-Refinement” discretization strategy is proposed for obtaining a mesh which is optimal for impact dynamics. In the case of oblique impact, the influence of Coulomb’s friction coefficient is investigated additionally. It shows that the numerical results are in good agreement with the experimental results for both normal and oblique impacts.
1. Introduction
For accurate modeling of multibody systems with impacts, there are two main ingredients that should be basically chosen while developing a contact formulation: one is the scheme to describe the kinematics and deformation of contact bodies and the other is the technique to enforce the nonpenetration constraint during contact.
The numerical approaches for kinematics and deformation description in contact/impact analysis can be divided into two categories: computational contact mechanics based on finite element method (FEM), as shown in, for example, [1–3], and approaches based on multibody system dynamics (MBS), as shown in, for example, [4–6]. In computational contact mechanics, nonlinear FEM is used to model the deformations, and usually the node-to-segment approach is used to discretize the contact surfaces. It is of major interest to predict the strain and stress field in the deformable bodies coming into contact. It is well suited for particularly high accuracy requirements yet causes very high computational effort. On the other hand, methods based on MBS provide a more efficient yet less accurate analysis of the contact/impact problems. These approaches emphasize capturing the effect of impact on the overall motion of the system for long simulation times. In order to gain computational efficiency, formulations based on the modal superposition and substructure methods have been used to reduce the number of degrees of freedom. Some simplifying assumptions are made in these methods: the contact interface is represented by two coincided points and a local reference frame rather than a time-varying area, or the contact force is modeled as a Hertz force element rather than distributed forces.
The other influential factor on the numerical result of dynamic simulation is the appropriate approach to describe impact. For contact/impact dynamics of discretized elastic body, two typical methods are presented to model contact: penalty method and Lagrangian method. A brief overview of these two formulations is presented as follows. In penalty method, the contact force is defined by a force function of local penetration at each contact point, as a spring-damper model. Various forms of force function are presented, developed, and used by a number of researchers [7–9], among others. This method has gained significant popularity because it does not bring extra dimensions to the dynamic equations. However, as the nonpenetration constraint is not precisely satisfied during the contact process using this method, the accuracy of numerical simulation depends on the choice of the coefficient of contact stiffness. The Lagrangian method for contact modeling, where nonpenetration condition is strictly and exactly satisfied and the contact force is represented by constraint matrix and Lagrange multipliers, has also attracted attentions from a lot of researchers [10–13]. Although the Lagrange multipliers increase the computation scale and bring numerical complexity, the constraints do not allow for any penetration at contact location, the Lagrangian method matches the physical reality, and not manually defined parameters are needed in the dynamic analysis.
For the experimental investigations of impact problems, the main difficulty is that the impact duration is very short and the impact responses have extremely high frequencies, so they are very difficult to measure. Early experimental investigation mainly focused on transient strain response and the measurement instruments are strain gauges [14]. Accelerometers were also used in some impact experiments to measure points located not close to the impacts [15, 16]. With the advancement of laser techniques, some impact experiments were conducted using Laser Doppler Vibrometers (LDVs) but mainly for one-dimensional rod impact or two-dimensional planar impact problems [17, 18].
In this work, a three-dimensional impact of a hemispherical-tip rod on a square plate, including normal impact and oblique impact, is investigated both numerically and experimentally. In the numerical approach, finite element method combined with the floating reference frame method is employed to describe the rigid-flexible coupling motion. To exclude the influence of several manually defined parameters and enforce the impenetrability constraint strictly, the Lagrangian method is applied to model the contacts. In the experimental study, strain gauges and LDVs are employed to measure the high-frequency impact responses. Great attention is given to investigate the influence of mesh size on the simulation accuracy, aiming to get a finite element mesh technique which ensures a convincing numerical solution. In our work, a “Prediction-Refinement” discretization strategy is proposed for obtaining a mesh which is optimal for impact dynamics, in the sense that the computational efforts involved are minimal under the accuracy requirement. In the case of oblique impact, the influence of Coulomb’s friction coefficient is investigated. The numerical results show good agreement with the experimental results for both normal and oblique impacts if the spatial discretization is conducted appropriately.
2. Contact/Impact Problem in Multibody Dynamics
In multibody dynamics with contacts/impacts, the rigid motion of flexible body is described by floating reference frame, and the deformation is modeled by finite element nodal coordinates or modal coordinates. This section provides the derivation of the kinematics and equations of motion of contact/impact problem in multibody system.
2.1. Kinematic Description
As shown in Figure 1, a deformable body Bi is discretized by lumped mass FEM; therefore, the mass of body is distributed to each finite element node. The inertial reference frame is represented by , and the floating reference frame attaching to the body is denoted by . The position vector of with respect to is denoted as . An arbitrary node position of flexible body is defined by (ignoring the body mark )where represents the position of node P in inertia frame, is the initial position, is the displacement in the floating frame, and is the transformation matrix from the floating frame to the inertial frame. In (1) and the following equations, the notation means quantity and is expressed in the floating frame.
Deriving (1) with respect to time yields the nodal absolute velocity, and deriving (1) twice results in the absolute acceleration:where represents the angular velocity of the floating frame, andEquations (2a) and (2b) can be expressed in matrix form as follows:where is the generalized velocity vector, and and are given byMatrix is a Boolean matrix related to node number, which leads to
2.2. Equations of Motion without Contact/Impact
Applying the principle of virtual power for the whole body, the dynamic equations of motion can be written as follows:where is the nodal external force vector and and are the damping matrix and stiffness matrix assembled by FEM.
Substituting (4a) and (4b) into (7) yieldswhere , , , and are given byFor a single flexible body, the coordinates are independent, and (8) can be written as follows:For the system containing bodies connected by kinematic joints, the equations of motion can be written as follows:where , , , and are matrix related to kinematic joint constraints, and represents the corresponding Lagrange multipliers.
2.3. Equations of Motion with Contact/Impact
There are several approaches to model the contact in the flexible multibody dynamics, including analytical solution, penalty method, and Lagrange multiplier method. Analytical solution such as St. Venant’s contact theory is limited to objects with simple geometry. In the penalty method, the contact force is defined by a force function of local penetration at the contact point, because the nonpenetration constraint is not precisely satisfied during the contact process; the accuracy of numerical simulation depends on the choice of the coefficient of contact stiffness. In contrast, the Lagrange multiplier method for contact modeling, where constraint equations are appended to dynamic equations to be solved together, reflects the nonpenetration condition without manually defined parameters. In this section, using Lagrange multiplier method, the contact constraint equations and equations of motion of system are derived.
In FEM, mostly the node-to-segment approach is used to discretize the contact surfaces, as shown in Figure 2. Possible contact points , of a pair satisfy the following conditions [10]:where and are the normal and tangential vector of contact surface and is denoted as the distance vector from to which is written as follows:
The position of point is interpolated by the surrounding nodes on the segment and is the interpolation function given bywhere and are the local coordinates of the segment.
When contact occurs, according to the nonpenetration condition in the contact point, the locking of the free motion in a normal direction is described by the following constraint equation on position level:Deriving (15) yields the constraint equation on velocity level:
According to (12a) and (12b), as is on the tangent plane, the second term in (16) disappears. Deriving (16) yields the constraint equation on acceleration level:
In the case of stick, the motion in the tangent direction is also locked. This results in the following constraints on velocity and acceleration level:
In (20) and (21), the second terms vanish according to (16), (18), and (19). Substituting (4a), (4b), and (13) into (17), (20), and (21), the normal and tangential constraints can be expressed aswhere represents the contact pair number, and are matrix related to normal contact constraints, and and are matrix related to tangential contact constraints.
If the contact pair is sliding, the friction force on body is modeled by Coulomb’s law:where is the normal contact force of this contact pair and is the Coulomb friction coefficient. Substituting (23) into (9a), (9b), (9c), and (9d), the generalized slip friction force can be obtained:
To describe the stick-slip phenomenon, a related transition law can be written as follows:
Assume, at a certain moment, there are contact pairs that are active, in which pairs are in stick state and contact pairs are in slip state. The equations of motion of the whole system with contacts and frictions can be assembled as follows:where describes the kinematic joint force, describes the normal contact force, and and describe the stick friction and slip friction, respectively. In each step, the contact state is determined by the contact detection and the stick or slip state is determined by the transition law written in (25a), (25b), (25c), and (25d). To solve (26), the explicit Runge-Kutta method is employed to accomplish the integration process.
3. Experimental Setup
A schematic diagram of the normal impact experiment is shown in Figure 3, and an overview of the experimental setup is shown in Figure 4. A cylindrical steel rod with hemispherical tip is used to strike an aluminum plate. The two colliding bodies are suspended by thin lines and are positioned horizontally by means of a spirit level. The rod just contacts the plate’s center and the rod’s axis is vertical to the plate’s surface when in equilibrium position. The plate keeps still until collision occurs. To determine the time of collision, we connect the two bodies to a direct current source, and the current signal comes out in the contact process. For the investigated duration of 5 ms, the motion of both bodies can be considered as a free horizontal motion.
The strains are measured with strain gauges. As shown in Figure 4, three gauges are bonded to the contact surface of the plate in three directions 15 mm away from the center point. The used signal amplifier is of type YE3817c made by Sinocera Piezotronics, Inc. In the measurement, the supplied constant voltage is set to 6 V and the amplifier gain is set to 500.
For the measurement of displacements and velocities, two Laser Doppler Vibrometers (LDVs) of type PSV-300F, made by Polytec GmbH, are used. The vibrometer utilizes an interferometric technique to measure vibrational signals. The measurement range of velocity is ±10 m/s, and the resolution reaches 10−6 m/s. The back point of the rod and the center point of the plate are measured in the normal impact experiment. The material and geometrical parameters of the two colliding bodies are tabulated in Table 1.
4. Comparisons between Measurements and Simulations
For simulating contacts using FEM, the modeling of contact itself as well as the consideration of the resulting wave propagation is closely related to the spatial discretization. Therefore, in order to obtain the convergence solution in FEM, the rod/plate case is simulated with different meshes and the simulation results are compared with the measurements.
For the main part of the two bodies, to get an accurate representation of wave propagation using the FEM, the smallest interesting wave length should be discretized by at least 20 nodes, as shown in, for example, [19]. Therefore, the maximum mesh size can be obtained:where is the highest frequency under consideration and is the lowest wave speed in the elastic materials. The longitudinal wave speed in the steel rod is
In the aluminum plate, the longitudinal and the lower transversal wave speeds can be written as follows:
To determine the value of the highest frequency , here, we perform a fast Fourier transform (FFT) of the velocity curve of the impact point P1, as shown in Figure 5. It can be seen that the frequency mainly concentrates in the range from 0 to 10 KHz, and the highest frequency of velocity response is up to 20 KHz. In this rod-plate case, the highest frequency is considered as 10 kHz firstly. Based on (27), the required maximal element size of the rod is 25 mm, and the plate is 15 mm.
For the discretization of the contact region, additional considerations are required. To sufficiently represent the deformation and stress distribution in the contact region, the contact region must be discretized in a much smaller size than the element length required for the wave propagation, as shown in Figure 6. To mesh the contact region, the maximum contact radius should be estimated beforehand. Here, we use an analytical method such as Hertz method to simulate the impact case firstly. It predicts a rough maximum force of about 800 N. And, from this, the contact radius a can be calculated by the Hertz law:where is the contact force, is the curvature radius of the contact surface, and, in this case, mm, .
As the contact radius is about 0.49 mm, to ensure that a certain amount of nodes are in contact, the element length in the contact region should be less than 0.49 mm. To investigate the influence of element size in the contact region, three meshes are used to discretize the bodies. In mesh 1, the contact element size is set to 0.4 mm amd mesh refinements are performed for meshes 2 and 3 in the contact region, while the elements in the main part keep unchanged. The parameters of the three meshes are listed in Table 2.
Figure 7(a) shows a comparison of the FEM simulation results with mesh 1 to mesh 3 for the velocity signal at the center point P1 of the plate. Small deviation is visible for the simulated response of mesh 1 compared with that of mesh 2 and the results of mesh 2 and mesh 3 are nearly the same. However, Figure 7(b) shows that the simulation using mesh 3 still does not agree well with the experiment although the element in the contact region is small enough. It means only refining the contact region still cannot lead to the accurate impact response. In this case, the impact process lasts about 1.7 ms according to the current signal. Obviously, the simulations using mesh 1 to mesh 3 have shorter impact duration and predict incorrect wave propagation. It means the element in the main part is not small enough to cover the high-frequency response caused by impact; thus, the mesh in the noncontact region has to be refined as well.
(a)
(b)
As shown in Table 2, mesh refinements are performed for mesh 4 and mesh 5 in the noncontact region, while the contact elements keep the same as mesh 2. The simulation results are plotted in Figure 8(a). Large differences exist among the results of mesh 3 to mesh 5. This phenomenon tells that the mesh size in the noncontact region also plays an important role in the simulation. For mesh 6, both the noncontact region and contact region are refined with smaller elements than mesh 5. Nearly no difference can be observed between simulations of mesh 5 and mesh 6. That is to say, the result of mesh 5 already reaches the convergence solution. In Figure 8(b), the result of mesh 5 is plotted together with the experiment, and the deviation nearly vanishes. The total node number of mesh 5 exceeds 90000 and a very small time step of is used to guarantee numerical stability, so the computational burden is very large. It takes over 20 hours to calculate the whole impact duration.
(a)
(b)
The other comparisons between simulation using mesh 5 and experiment, including velocity of point P2 and strain of point S2, are shown in Figure 9. It shows that simulation using mesh 5 is consistent with the measurement.
(a)
(b)
In this section, the rod-plate impact problem is simulated using 6 meshes and finally the numerical result agrees well with the experimental result. It means that the finite element mesh is influential to the simulation result. Here, a “Prediction-Refinement” mesh strategy is shown in Figure 10, aiming to obtain an optimal mesh for impact dynamics. When an impact problem is investigated, the following steps for spatial discretization can be considered.
(1) Prediction. For discretization of the contact region, an analytical method such as Hertz model is used to estimate the maximum contact radius, and the element length in the contact region should be less than the contact radius to ensure that a certain amount of nodes are in contact. And the discretization of the noncontact region is calculated from the lowest wave speed (material) and the highest domain frequency under consideration (measurement or experience).
(2) Refinement. The predicted discretization still cannot lead to an accurate result. The mesh refinement should be performed both in the contact region and in the noncontact region until the result reaches convergence. In this way, the optimal mesh is obtained, in the sense that the computational efforts involved are minimal under the accuracy requirement.
5. Oblique Impact
To investigate the tangential motion during the contact, an oblique impact experiment is performed. The experimental equipment is the same as in the normal case. The vertical view of the experimental setup is shown in Figure 11. The axis of the rod points to the plate center and the angle between the rod and the plate is 45°. To reflect the role of the friction, the velocity of the point on the plate side is measured by the LDV. The measurement procedure for the oblique impact is similar to the procedure for the normal impact.
In the simulation, the mesh size used here is the same as the mesh 5 in the normal case. The friction model is implemented to represent the tangential force, and the friction coefficient is set as 0.35 according to measurement. In order to investigate the influence of the friction coefficient on the tangential motion, the friction coefficients of 0.15 and 0.55 are also simulated.
Figure 12 shows the computed and measured velocities of the plate and the rod. Even though the tangential force is smaller than the normal contact force and acts only during a small time period, it has a measurable effect on the motion of the contact bodies. Different friction coefficients in the simulation lead to dissimilar numerical results. The tangential motion of the plate is directly influenced by the friction force. The behavior of the rod after the oblique impact is a combination of lateral and rotational rigid body motion superposed by the elastic vibrations. The friction force influences both the lateral tangential velocity and the angular velocity; thus, the simulated velocity of the measured point on the rod varies due to different friction coefficients.
(a)
(b)
In Figure 12(a), a chattering phenomenon happens in the experimental curve, and it is mainly because the location of the measured point changes due to the rigid motion caused by the impact. The numerical results agree well with the experimental results when the friction coefficient is set properly.
6. Conclusion
In this paper, the impacts between a steel rod and an aluminum plate are studied numerically and experimentally. In the numerical approach, FEM is used to discretize the contact bodies to describe the deformation precisely combined with the floating reference frame method to describe the rigid motion. The Lagrangian method is applied to model the normal contacts to enforce the impenetrability constraints strictly and exclude the influence of some manually defined parameters such as contact stiffness in penalty method. To model the tangential friction, the tangential contact is also modeled using Lagrangian method and the sliding friction is described by Coulomb’s friction law. In the experimental study, strain gauges and LDVs are employed to measure the high-frequency impact responses.
Firstly, the normal impact case is investigated, and detailed comparative studies between the numerical and experimental results show that the simulation accuracy is greatly influenced by the spatial discretization. This work presents a “Prediction-Refinement” discretization procedure to obtain a mesh which is optimal for impact dynamics using FEM, in the sense that the computational efforts involved are minimal under the accuracy requirement. For the discretization of the contact region, the contact radius can be predicted by an analytical method such as Hertz model in advance. The expected contact region has to be discretized by sufficient elements to ensure that a number of nodes are in contact in order to describe the contact accurately. The discretization of the main part is also essential to the numerical results. Because the impact leads to very high-frequency vibration in the whole body, the required mesh size for correct computation of the wave propagation can be calculated from the lowest wave speed and the highest considered frequency. Then the mesh should be refined both in the contact region and in the noncontact region until the numerical result reaches convergence. The simulation coincides with the measurement only when the contact bodies are discretized appropriately.
In the case of oblique impact, the experiment shows that friction has a measurable influence on the tangential motion. The value of friction coefficient greatly affects the resultant responses. Good agreements between the simulations and measurements can be achieved if the friction coefficient is chosen reasonably.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research work is supported by the National Science Foundation of China (no. 11132007 and no. 11202126), for which authors are grateful.