Multibody Syst Dyn(2010)24:103–122
DOI10.1007/s11044-010-9209-8
On the contact detection for contact-impact analysis
in multibody systems
Paulo Flores·Jorge Ambrósio
Received:28May2008/Accepted:2May2010/Published online:20May2010
©Springer Science+Business Media B.V.2010
Abstract One of the most important and complex parts of the simulation of multibody sys-tems with contact-impact involves the detection of the preci instant of impact.In general, the periods of contact are very small and,therefore,the lection of the time step for the integration of the time derivatives of the state variables plays a crucial role in the dynamics of multibody systems.The conrvative approach is to u very small time steps through-out the analysis.However,this solution is not efficient from the computational view point. When variable time-step integration algorithms are ud and the pre
impact dynamics does not involve high-frequencies,the integration algorithms may u larger time steps and the contact between two surfaces may start with initial penetrations that are artificially high. This fact leads either to a stall of the integration algorithm or to contact forces that are physically impossible which,in turn,lead to post-impact dynamics that is unrelated to the physical problem.The main purpo of this work is to prent a general and comprehensive approach to automatically adjust the time step,in variable time-step integration algorithms, in the vicinity of contact of multibody systems.The propod methodology ensures that for any impact in a multibody system the time step of the integration is such that any initial penetration is below any prescribed threshold.In the ca of the start of contact,and after a time step is complete,the numerical error control of the lected integration algorithm is forced to handle the physical criteria to accept/reject time steps in equal terms with the numerical error control that it normally us.The main features of this approach are the simplicity of its computational implementation,its good computational efficiency,and its ability to deal with the transitions between non-contact and contact situations in multibody dynamics.A demonstration ca provides the results that support the discussion and show the validity of the propod methodology.
P.Flores( )
Departamento de Engenharia Mecânica,Universidade do Minho,Campus de Azurém,4800-058
Guimarães,Portugal
e-mail:pflores@dem.uminho.pt
J.Ambrósio
Departamento de Engenharia Mecânica,Instituto Superior Técnico,IST/IDMEC,Av.Rovisco Pais,1, 1049-001Lisbon,Portugal
e-mail:jorge@dem.ist.utl.pt
104P.Flores,J.Ambrósio Keywords Contact detection·Contact-impact analysis·Time integrators·Integration error control·Variable time step·Multibody dynamics
1Introduction
The classical problem of the contact mechanics is still an open issue in engineering appli-cations.In particular,the contact-impact modeling and analysis in multibody dynamics has received a great deal
of attention over the past few decades and still remains an activefield of rearch and development[1–6].Contact events happen frequently in multibody systems and in many cas the function of mechanical systems is bad on them[7–15].In general, the motion characteristics of a multibody system are significantly affected by contact-impact phenomena.Impact is a complex physical phenomenon for which the main characteristics are a very short duration,high force levels,rapid energy dissipation,and large changes in the velocities of the bodies[16].Inherently,contact implies a continuous process which takes place over afinite time.Other effects directly related to the impact phenomena are tho of vibration propagation through the system,local elastic/plastic deformations at the contact zone,and frictional energy dissipation[17–26].Impact is a prominent phenomenon in many mechanical systems such as mechanisms with intermittent motion and mechanisms with clearance joints[27–31].In a broad n,the contact-impact modeling in multibody systems consists of two major steps,namely,the contact detection and contact respon.
The subject of development of contact detection problem is a quite challenging and actual problem in variousfields such as,discrete element methods[32],robotics[33],or vehicle systems[34].From the modeling methodology point of view,veral different methods have been developed.Carsten and Wr
iggers[35]prented an explicit multibody contact algo-rithm where the contact detection issue was also studied using a predictor-corrector scheme. An iterative form of the propod scheme was also ud to reduce the computational effort. One of the most robust and well-known methods for contact detection of complexly shaped bodies was propod by Hippmann[36].This algorithm,referred to as polygonal contact model,is bad on reprentation of the body surfaces by polygon meshes and the contact force evaluation is done using an elastic foundation model.This approach has been ud by other rearchers[37,38].He et al.[39]prented a multigrid contact detection method, where the multigrid idea was integrated with contact detection problems.Wellmann et al.
[40]developed an efficient contact a contact detection algorithm for super-ellipsoids bad on the common-normal concept.The problem of contact detection is formulated as2D un-constrained optimization problem that is solved by a combination of Newton’s method and a Levenberg–Marquardt method.More recently,Studer et al.[41]extended the modern time-stepping algorithms to include a step-size adjustment and extrapolation for contact detection in nonsmooth dynamics.Portal et al.[42]prented a methodology for contact detection between convex quadric surfaces using its implicit equations.This methodology was im-plemented in a multibody dynamics code in order to simulate the interpenetration between mechanical systems,particularly,the simulation of collisions wit
h automotive vehicles and other road urs,such as cars,motorcycles,and pedestrians.The contact detection of two bodies was formulated as a convex nonlinear constrained optimization problem that is solved using two methods,an Interior Point method(IP)and a Sequential Quadratic Programming method(SQP).
From the modeling methodology point of view,veral different methods have been in-troduced to model the contact respon in multibody systems.As a rough classification,they may be divided into contact force bad methods[43]and methods bad on geometrical
On the contact detection for contact-impact analysis105 constraints[44],each of them showing advantages and disadvantages for each particular application.
Contact force approaches,commonly referred as penalty or compliant methods,own their importance in the context of multibody systems with contacts to their computational sim-plicity and efficiency[17].In the methods,the contact force is expresd as a continuous function of penetration between contacting bodies.One of the main drawbacks associated with the force models is the difficulty to choo contact parameters such as the equiva-lent stiffness or the degree of nonlinearity of the penetration,especially for complex contact scenarios and nonmetallic materials[
45].The penalty formulations can be understood as if each contact region of the contacting bodies is covered with some spring-damper elements scattered over its surfaces.The normal force,including elastic and damping,prevents ,no explicit kinematic constraint is considered but simply force reaction terms are ud.The magnitudes of stiffness and deflection of the spring-damper elements are com-puted bad on the penetration,material properties and surface geometries of the colliding bodies.In the work by Khulief and Shabana[28,29]the required parameters for reprenting contact force laws are obtained bad on the energy balance during contact.This formulation us a force-displacement law that involves determination of material stiffness and damping coefficients.In the work by Lankarani and Nikravesh,[46]two continuous contact force models are prented for which unknown parameters are evaluated analytically.In thefirst model,internal damping of bodies reprents the energy dissipation at low impact velocities. However,in the cond model,local plasticity of the surfaces in contact becomes the dom-inant source of energy dissipation.Dias and Pereira[47]described the contact law using a continuous force model bad on the Hertz contact law with hysteresis damping.The effect and importance of structural damping schemes inflexible bodies were also considered.Hunt and Crossley[48]obtained a model for computing the stiffness coefficient from the energy balance relations.Bad on the Hunt and Crossley approach Lankarani and Nikravesh[43] further extend the contact model with hysteresi
s damping.In their approach,the damping force is a linear function of the elastic penetration which is estimated from the energy dissi-pated during impact.The effect of friction in this approach is often taken into consideration by using a regularized Coulomb friction model.An overview of different models of friction together with fundamentals can be found in Oden and Martins[49]and Feeny et al.[50].
The complementarity formulations associated with the Moreau’s time-stepping algorithm for contact modeling in multibody systems have ud by many rearchers[3,10,32,38, 41].Assuming that the contacting bodies are truly rigid,as oppod to locally deformable or penetrable bodies as in the penalty approaches,the complementarity formulations resolve the contact dynamics problem by using the unilateral constraints to compute contact im-puls or forces to prevent penetration from occurring.Thus,at the core of the complemen-tarity approach is an explicit formulation of the unilateral constraints between the contacting rigid bodies[51].The basic idea of complementarity in unilateral multibody systems can be stated as for a unilateral contact either relative kinematics is zero and the corresponding constraint forces are zero or vice versa.The product of the two groups of quantities is al-ways zero.This leads to a complementarity problem and constitutes a rule which allows the treatment of MBS with unilateral constraints[52–55].One of thefirst published works on the
complementarity problems is due to Signorini[56],who introduced an impenetrability condition in the form of a linear complementary problem.Later,Moreau[57]and Pana-giotopoulos[58]also applied the concept of complementarity to study nonsmooth dynamic systems.Pfeiffer and Glocker[7]extended the developments of Moreau and Panagiotopou-los to multibody dynamics with unilateral contacts,being the complementarity considered of paramount importance.Indeed,complementarity problems proved to be a very uful way to formulating problems involving discontinuities[59–62].
106P.Flores,J.Ambrósio In a dynamic simulation,it is very important tofind the preci instant of transition be-tween the different states,that is,the transition between contact and noncontact situations. Especially when continuous contact force models are ud,such as the one propod by Lankarani and Nikravesh[43],if the instant of the start of contact is not detected properly the initial contact force may become abnormally large due to the unphysical high initial pen-etrations between the impacting surfaces.This numerical abnormality leads to an artificial increa of the system energy and,eventually,to the stall of the integration process,when variable time-step integration algorithms are ud.The avoidance of this problem requires a clo monitoring of the numerical procedure to continuously detect and analyze all situ-ations.Otherwi,the errors may buildup and thefinal results are meaningless.Thus,the main purpo of this work is to prent a general and com
prehensive methodology to deal with the detection of the preci instant of contact in multibody dynamics and to propo actions at the level of the integration algorithm that,without interfering with its mathemati-cal structure,allows controlling time steps bad on physical reasoning as a complement of the time-step control inherent to all variable time-stepping integration algorithms.
2Model for contact forces
In order to evaluate efficiently,the contact-impact forces resulting from collisions in multi-body systems special attention must be given to the numerical description of the contact force model.Information on the impact velocity,material properties of the colliding bod-ies,and geometry characteristics of the contact surfaces must be included into the contact force model.Due to its simplicity and ability to characterize the contact phenomena the contact forces are reprented,in this work,using a continuous force model bad on a penalty formulation[43].Further,it is important that the contact force model can add to the stable integration of the multibody system equation of motion.This contact force model is bad on the Hertz elastic contact law,being the hysteresis damping function incorporated to reprent the energy dissipated during the impact.Lankarani and Nikravesh[43]suggest parating the normal contact force into elastic and dissipative components as
f N=Kδn+D˙δ(1) where K is the generalized stiffness constant andδis the relative normal indentation be-tween the bodies.The exponent n is t to1.5for the cas where there is a parabolic dis-tribution of contact stress,as in the original work by Hertz[63].Convenient expressions for the contact force bad on experimental or numerical work u n=1.5,for metallic materials,and other exponents for other materials such as glass or polymers.Although such penalty formulations have the same form of(1)and they po the same numerical chal-lenges if the exponent n=1.5or the distribution of the contact stress is not parabolic, they must not be confud with the Hertz theory.The generalized parameter K is dependent on the material properties and the shape of the contact surfaces.For two spheres in contact, the generalized stiffness coefficient is function of the radii of the spheres i and j and the material properties as[64],
K=
4
3(σi+σj)
R i R j
R i+R j
1
2
(2)
On the contact detection for contact-impact analysis107 where the material parametersσi andσj are given by
σk=1−ν2
k
E k
,(k=i,j)(3)
and the quantitiesνk and E k are the Poisson’s ratio and the Young’s modulus associated with each sphere,respectively.
In(1),the quantity D is a hysteresis coefficient and˙δis the relative normal impact ve-locity.The hyster
esis coefficient is written as a function of penetration as
D=χδn(4) in which the hysteresis factorχis given by
χ=3K(1−c2
e
)
4˙δ(−)
(5)
being˙δ(−)the initial impact velocity.By substituting(5)into(4)and the result into(1),the normal contact force isfinally expresd as
f N=Kδn
1+
3(1−c2
e
)
4
˙δ
˙δ(−)
(6)
where the generalized parameter K is evaluated by(2)for sphere to sphere contact,or by
similar expressions for the contact of other types of geometry and c e is the restitution coef-
ficient.Therefore,it is crucial,for the correct u of the continuous force model,the exact
identification of the initial contact velocity˙δ(−)and to start the analysis of the contact pe-riod with a null penetrationδ,which in numerical terms means a penetration smaller than a
predefined ,δ<ε.
3Contact detection methodology
When a system consists of fast and slow moving components,that is,the eigenvalues are
widely spread,the system is designated as being stiff[65].Stiffness in the system equations
of motion aris when the gross motion of the overall multibody system is combined with the
nonlinear contact forces that lead to rapid changes in velocity and accelerations.In addition,
when the equations of motion are described by a coupled t of differential and algebraic
equations,the error of the respon system is particularly nsitive to constraints violation,
which inevitably leads to artificial and undesired changes in the energy of the system.Yet,by
applying a stabilization technique,the constraint violation can be kept under control[66].
During the numerical integration procedure,both the order and the step size of the inte-
gration algorithms are adjusted to keep the error tolerance under control.In particular,the
variable step size of the integration scheme is a desirable feature when integrating systems
that exhibit different time scales,such as in multibody systems with impacting bodies[67].
Thus,large steps are generally taken when the motion of the system does not include contact
forces but,when impact occurs,the step size is decread substantially due to the inclusion
of high frequency contents in the system respon and not becau the amount of penetration
obrved between two contact surfaces is larger or smaller.
One of the most critical aspects in the dynamic simulation of the multibody systems with
collisions is the detection of the preci instant of contact.In addition,the numerical model