A Differential Algebraic Method to Approximate Nonsmooth Mechanical Systems by Ordinary Differential Equations

and Applied Analysis Hindawi Publishing Corporation http://www.hindawi.com Volume 2013 ISRN Applied Mathematics Hindawi Publishing Corporation http://www.hindawi.com Volume 2013 Hindawi Publishing Corporation http://www.hindawi.com Volume 2013 International Journal of Combinatorics Hindawi Publishing Corporation http://www.hindawi.com Volume 2013 Journal of Function Spaces and Applications International Journal of Mathematics and Mathematical Sciences Hindawi Publishing Corporation http://www.hindawi.com Volume 2013


Introduction
Mechanical systems involving dry friction and rigid unilateral contact are usually described as differential inclusions (DIs).Conventional approaches for simulating those nonsmooth systems can be broadly categorized into two types: regularization approaches and hard-constraint approaches [1,2].
In regularization approaches, also referred to as penaltybased approaches [3,4], the discontinuous force laws of dry friction and rigid unilateral contact are approximated by using continuous functions.For example, some previous friction models [5][6][7][8][9] and contact models [10][11][12][13] can be viewed as approximations of dry friction and rigid unilateral contact, respectively.Physical meanings of such approximations can be usually interpreted as relaxation of constraints, that is, compliance that replaces rigid constraints between force and motion.By employing those models, the equations of motion of nonsmooth systems can be written by ODEs.As a result, discontinuous natures of the systems are lost, and consequently, some unrealistic behaviors can be produced.For examples, Dahl model [5] and LuGre model [6] can produce positional drift in the static friction state.
In hard-constraint approaches, rigid bodies are considered strictly impenetrable to each other.One major way of this approach is to discretize the equation of motion by backward Euler-like methods.The discretized equation is regarded as an algebraic equation, which is then solved numerically [14][15][16][17][18][19][20][21][22] or analytically ([23, Section III.A], [24,Section 1.4.3.2], and [25]).Another type of approach (e.g., [26]) is to describe a system as an ODE in every period between discontinuous events such as transitions between static and kinetic friction states.It is easy to see that such a scheme is not suitable when too many discontinuous events occur.
This paper introduces a new method to approximately describe nonsmooth mechanical systems by ODEs.This method is derived based on the observations that DIs describing nonsmooth mechanical systems can be approximated by differential algebraic inclusions (DAIs) and that those DAIs are equivalently rewritten as ODEs.In contrast to conventional regularization methods, this method preserves the intrinsic nature of discontinuity in those systems.This method is illustrated by some examples, of which simulation results are obtained through the fourth-order Runge-Kutta (RK4) method.
The rest of this paper is organized as follows.Section 2 gives some mathematical preliminaries to be used in subsequent sections.Section 3 overviews previous approximation methods for dry friction and rigid unilateral contact.Section 4 gives the main contribution of the work.Section 5 provides two example applications of the new method.Finally, concluding remarks are given in Section 6.

Mathematical Preliminaries
For the discussion throughout this paper, this section introduces three functions: sgn, sat, and dio.Some theorems regarding those functions are also presented.In the rest of this paper, R denotes the set of all real numbers and R + denotes the set of all nonnegative real numbers.The symbol 0 denotes the zero vector of an appropriate dimension.First, let us define the signum function sgn : R  → R  and the unit saturation function sat : R + × R  → R  as follows: sat where  ∈ R  ,  ∈ R + , and ‖⋅‖ denotes the vector two-norm.
If  = 1, the sgn() and sat(, ) can be depicted as Figures 1(a) and 1(b), respectively.The following theorem is useful to rewrite the DIs involving sgn as ODEs involving sat.
Theorem 1.For ,  ∈ R  and  ∈ R + , the following relation holds true [27,28]: Next, let us define the "diode" function dio : R + → R + as follows: where  ∈ R + .The following theorem is useful to rewrite DIs involving dio by ODEs.
Theorem 2. For  ∈ R and  ∈ R + , the following relation holds true: Proof.A proof can be given as follows: The graphs of dio() and max(0, −) are illustrated as Figures 1(c) and 1(d), respectively.
It must be noted that Theorems 1 and 2 are special cases of the following relation, which has been used in, for example, [24,Appendix A.3], [29, equation (2)], and [30, equation (4)]: Here,  ∈ R  ,  ∈  ⊂ R  ,  is a closed convex set,   () is the normal cone to the set  at , and prox(, ) is the "proximal point" function defined as follows: Theorems 1 and 2 can be obtained by using the relation (8) with  = { ∈ R  | ‖‖ ≤ } and  = R + , respectively.

Previous Approaches
3.1.Dry Friction.Let us consider the situation where a rigid mass  > 0, of which the position is  ∈ R  ( ∈ {1, 2}), is sliding on a fixed surface.Let us assume that it is subject to the dry friction force  ∈ R  and an external force   ∈ R  .Then, the equation of motion of the mass is described as the following DI: where and  > 0 is the magnitude of kinetic friction force.
(Common definitions of dry friction assume that the static friction force can be larger than the kinetic friction force.This paper leaves this out of consideration and assumes that the maximum static friction force is equal to the kinetic friction force.)The direct integration of ( 10) and ( 11) is difficult since the value of sgn( ṗ ) is not determined at ṗ = 0, according to the definition (1) of sgn.Some previous friction models can be viewed as approximations of (11).One simple way is to employ a threshold velocity [8,31] below which the velocity is considered zero.This method may be useful to avoid the discontinuity in (11), but the nonphysical threshold can produce unrealistic artifacts.Another way is to employ a new state variable which usually can be interpreted as the displacement of a viscoelastic element.For example, LuGre friction model [6] without Stribeck effect can be described as follows: where  ∈ R  is the new state variable,  > 0 is a sufficiently large constant, and ,  > 0 are constants appropriately chosen to suppress the oscillation in .Dahl friction model [5] is a special case of LuGre friction model with  =  = 0.
A disadvantage of those two models is that they produce unbounded positional drift in the static friction state under oscillatory external force even smaller than the maximum static friction force [23,32].Other types of regularized friction models are proposed by Kikuuwe et al. [23,Section III.C] and Bastien and Lamarque [33] based on Backward-Euler method and by Kikuuwe and Yamamoto [34] based on a modified Runge-Kutta method.A downside of their models is that they restrict the choice of methods for time integration.
In hard-constraint approaches, the equations of motion are discretized along time by Euler-like methods.Those discretized equations are usually formulated into complementarity problems, which are then numerically solved.The literature includes some complementarity formulations of dry friction in one-dimensional space [16,22] and in multidimensional space [15,[17][18][19].One exception is Kikuuwe et al. 's approach [23, Section III.A], in which the discretized equation in a very simple case is analytically solved by the application of Theorem 1 in the present paper.

Rigid Unilateral Contact.
Let us consider that the onedimensional system is composed of a rigid mass , of which the position is  ∈ R, and a fixed rigid wall whose position coincides with the origin.The rigid mass is subject to an external force   ∈ R.Then, the equation of motion of the rigid mass is described as the following DI: where The integration of ( 13) and ( 14) is also difficult due to dio(), whose value is not determined at  = 0.
One of the trivial methods to approximately realize the contact force  in ( 14) is as follows [24,34,35]: where  is a sufficiently large positive constant and  is a positive constant large enough to dampen the oscillation in .This force law can be viewed as a linear viscoelastic contact model with the stiffness  and the viscosity .As pointed out in [13,25], one drawback of ( 15) is that it produces an unnatural sucking force toward the wall.This drawback may be overcome by using the following slightly different one: However, both (15) and ( 16) are discontinuous with respect to  and ṗ .Thus, they are not suitable for the use with common ODE solvers.
As another example, the nonlinear viscoelastic contact model proposed by Hunt and Crossley [13] can also be viewed as an approximation of rigid unilateral contact.This model was extended in [11,12,36] and empirically validated in [10,37,38].This model is continuous with respect to  and ṗ , but it can also produce unnatural sucking force when ṗ is large.
In hard-constraint approaches for rigid unilateral contact, the equations of motion are usually discretized by Euler-like methods and then solved numerically [14,16,17,22,39].A different approach is in [24, Section 1.4.3.2] and [25] where the discretized equations in very simple cases are solved analytically.Those methods can be used only with Euler-like methods.

New Method
In this section, new ODE approximations are introduced for ((10), (11)) and (( 13), ( 14)).Based on those simple approximations, a general procedure is presented for approximating nonsmooth mechanical systems involving many rigidunilateral and dry-frictional contacts.

Dry Friction.
The new approach for approximating (11) is motivated by Kikuuwe et al. 's work [23].Their work (specifically, model-C in [23]) provides an idea to approximate (11) by the following DAI: Here,  ∈ R  is a state variable newly introduced,  > 0 is a sufficiently large constant, and  > 0 is a constant appropriately chosen to suppress the oscillation in .A physical interpretation of the approximation ((17a) and (17b)) can be illustrated as Figure 2. A friction force described by  sgn( ṗ − ȧ ) acts on a massless object whose velocity is ṗ − ȧ , and a viscoelastic element with the stiffness  and the viscosity  produces the force  in (17b), which exactly balances the friction force.

Massless object 𝐾 𝐾𝛽 𝑓
In Kikuuwe et al. 's method, (17a) is discretized by Backward-Euler method; for example, ȧ is replaced by (  −  −1 )/, where  denotes the timestep size and the subscripts denote time indices, and then it is analytically solved with respect to   by using Theorem 1.In Bastien and Lamarque's model [33], a set of inclusions and equations with similar form to ((17a) and (17b)) are also discretized by Backward-Euler method and then analytically solved.
The observation that motivated the new approach is that ((17a) and (17b)) can be solved without using the Backward-Euler method.By the direct application of Theorem 1, ((17a) and (17b)) can be equivalently rewritten as the following ODE: As far as the authors are aware, the literature includes no computational methods making use of the equivalence between DAIs of the form of ((17a) and (17b)) and ODEs of the form of ((18a) and (18b)).
After replacing ( 11) by (18b) and appending (18a) to the state-space model, the system (10) and ( 11) is approximated by the following ODE: Figure 3 shows the simulation result by using the ODE (19) with RK4.To illustrate the advantage of this method, it also presents the result of LuGre model ((12a) and (12b)) combined with (10).In the simulation, an external force   was chosen as which is, after  = 2 s, oscillatory below the static friction level  = 0.5 N. As shown in Figure 3(b), LuGre model produces unrealistic positional drift, which has been known in the literature (e.g., [23,32]), while the presented method (19) does not.This implies that ((18a) and (18b)) is a better approximation of ( 11) than ((12a) and (12b)).
It should be mentioned that ((18a) and (18b)) is derived by relaxing the rigid constraint between  and ṗ in (11) by introducing an auxiliary variable  that has its own dynamics.In this sense, the proposed method may be viewed to be similar to Baumgart's method [40], in which constraints are relaxed to improve the numerical stability of the solutions of ODEs.One of the concerns about DIs is the existence and uniqueness of their solutions, as discussed by Bastien and Lamarque [33].As for the case of ((17a) and (17b)), on the other hand, it is clear because ((17a) and (17b)) is equivalent to the ODE ((18a) and (18b)).

Rigid Unilateral Contact.
The new approach for approximating ( 14) is a modification of the work by Kikuuwe and Fujimoto [25].In their approach, ( 14) is approximated by the following DAI: where  and  are appropriate positive constants and  ∈ R is a state variable newly introduced.A physical interpretation of ((21a) and (21b)) is illustrated as Figure 4. Here, a massless object whose position is  +  is connected to the mass through a viscoelastic element with the stiffness  and the viscosity .Due to the contact, the contact force dio( + ) acts on the massless object and it balances the force  from the viscoelastic element.In Kikuuwe and Fujimoto's work, ((21a) and (21b)) was discretized by Backward-Euler method and then analytically solved by the application of Theorem 2. Unfortunately, ((21a) and (21b)) cannot be rewritten into an ODE because ė cannot be obtained explicitly.
The new approach presented here is to add another term  ė to the argument of dio(⋅), which yields the following DAI: where  > 0 is another appropriate constant.By using Theorem 2, (( 22) and ( 23)) can be equivalently rewritten as the following ODE: The equivalence between DAIs of the form of (( 22) and ( 23)) and ODEs of the form of (( 24) and ( 25)) has not been pointed out in the literature either.One can see that (( 24) and ( 25)) is continuous with respect to , ṗ , and , and it does not produce sucking force because the right-hand side of ( 25) is always positive.This features in contrast to the conventional methods ( 15) and ( 16), which are discontinuous, and to Hunt and Crossley's model [11][12][13], which produces a sucking force.It is also easy to see that (( 22) and ( 23)), or equivalently (( 24) and ( 25)), has a unique solution.
One possible interpretation of (( 22) and ( 23)) and its equivalent expression (( 24) and ( 25)) can be explained by defining ẽ ≜ + ė . By using ẽ, (( 22) and ( 23)) can be rewritten as follows: where L denotes the Laplace transform.By noting the similarity between ((26a) and (26b)) and ((21a) and (21b)), one can see that force  in (26b) can be interpreted as a lowpass filtered viscoelastic force although it does not exist in the real world.When  = , ((26a) and (26b)) is equivalent to ((21a) and (21b)) with  = 0, which produces a perfectly elastic force.To preserve the effect of the viscous force, it is presumable that  should be set smaller than , although any tuning guidelines are not obtained yet.By replacing ( 14) by (25) and appending (24) to the statespace model, the system ( 13) and ( 14) is approximated by the following ODE: A set of numerical simulation of the ODE ( 27) was performed with different  values and a fixed  value.Figure 5 shows that the bouncing motion becomes smaller as  decreases.This is consistent with the interpretation based on ((26a) and (26b)), which implies that a smaller  strengthens the viscous effect in a high-frequency region.Detailed analysis on the relation between the parameter values and the achieved coefficient of restitution is left outside the scope of this paper.What can be said is that the coefficient of restitution can be adjusted by appropriate choices of  and  on a trial-and-error basis.
where  = [   ,   ]  ,   ∈ R 2 , and  is the friction coefficient between the mass and the surface.
It must be noticed that (28) includes a multiplication of dio(⋅) and sgn(⋅).To approximate this, one must replace dio(⋅) first and then replace sgn(⋅) because the replacement of sgn(⋅) involves its multiplicative factor ( in (11) ) while that of dio(⋅) can be done independently.In conclusion, the DI ( 28) can be approximated by the following ODE: where and the parameters ,  1 ,

General
Procedure.Now we are in a position to present the main contribution of the work.A mechanical system can be generally described by a DI in the following form: where  ∈ R  is the state vector of the system.Here, Φ is a function that contains dio(⋅) and sgn(⋅) in several places and may also contain single valued functions.With this procedure, the nonsmooth system (31) is approximated by the following ODE: where Φ(,  1 , . . .,   ,  1 , . . .,   ) is the function Φ() in which the aforementioned replacements are made.
The presented procedure cannot apply if the function Φ() includes a sgn(⋅) whose multiplicative factor involves discontinuous functions other than dio(⋅) and if   () are not guaranteed to be nonnegative.The authors, however, are not aware of nonsmooth mechanical systems that must be described by such Φ() functions.

Example I: A Rolling Sphere with Collision and Slip.
The presented approach is now illustrated by an example problem.Let us consider a system in which a spherical object with a uniform mass density falls onto a fixed rigid surface, as shown in Figure 6.The surface includes the origin and is perpendicular to the -axis.This example is also introduced in [34] and a similar example is employed in [11].Let  ∈ R 3 be the position of the gravity center of the object,  be the unit quaternion representing the attitude of the object, and  ∈ R 3 be the angular velocity of the object.Let  and  be the radius and the mass of the object, respectively, and let  be the friction coefficient between the object and the surface.
Then, the equations of motion of the object can be described as the following DI:
According to the procedure presented in Section 4.4, the DI (33) can be approximated by an ODE in the following procedure.First, one should replace dio(  − ) by where  ∈ R and where  ∈ R 2 and  2 and  2 are positive constants appropriately chosen.Finally, ODEs defining the behaviors of  and  should be appended to (33).Then, ( 33) is approximated by the following ODE: Figure 7 shows the result of the simulation by using (37) with RK4.Here,  1 and  2 are set as high as possible to achieve small penetrations during collisions, and  1 ,  2 , and  are chosen based on some trials and errors.Figures 7(a In Figure 7(e), one can find small penetrations produced by the approximation.Moreover, in Figure 7(f), one can see impulse oscillations after collisions, which are also consequences of the approximation.Despite these small artifacts, the overall shapes of the graphs in Figures 7(a) to 7(d) are close to the expected behaviors of the original DI (33).

Example II: Multiple Frictional-Unilateral Contacts.
Next example is the application of the presented method to a system involving many frictional contacts interacting with one another.Let us consider a planar system illustrated in Figure 8, which consists of a conveyor moving at a constant velocity , a spring with the stiffness   , two rigid objects  1 and  2 , and a rigid vertical wall.The object  1 can move freely in the horizontal direction and is subject to the elastic force from a spring   in the vertical direction.It is assumed that the objects do not rotate.The coefficients of friction between the wall and  1 , between  1 and  2 , and between  2 and the conveyor are  1 ,  2 , and  3 , respectively.The state vector is defined as denote the positions of  1 and  2 , respectively.The object  1 is regarded as being at [0, 0]  when it is in contact with the wall and the spring balances the gravity.The object  2 is regarded as being at [0, 0]  when it is in contact with the conveyor and the object  1 being at its [0, 0]  .Then, the state-space model of the system can be described as the following DI:

dio (𝑝
where and   and   ( ∈ {4, 5, 6}) are appropriate constants.Finally, the differential equations defining the behaviors of   and   ( ∈ {1, 2, 3}) should be appended to (38).Then, (38) is approximated by the following ODE: A numerical simulation was performed by using the ODE (45) with RK4.The results are shown in Figure 9. Here, again,   are set as high as possible to achieve small penetrations during collisions, and   and   are chosen based on some trials and errors.The time periods indicated by the gray regions are those in which the objects  1 and  2 are in contact to each other.Figures 9(a) and 9(b) show the horizontal bouncing motion of  1 and  2 , which eventually converges.Figure 9(c) shows the vertical motion of  1 , which exhibits nonsmooth changes in the velocity during the contact with  2 , being influenced by the friction force.The vertical position of  1 does not converge to zero because of the static friction forces from the wall and  2 .Figure 9(d) shows the vertical motion of  2 , which determines the normal force from the conveyor to  2 .It properly shows the influence on the normal force from the friction force acting on the side face.
Also in this simulation, one can observe small penetrations at the time of collisions in Figures 9(a) and 9(d).In addition, in Figure 9(b), one can find some impulsive responses in ṗ 1 , which are also caused by the approximation.Despite these artifacts, one can see that the approximation (45) appropriately simulates the overall behavior of the original DI (38).

Example III: Periodic Motion.
This section shows the application of the proposed method to a system exhibiting periodic motion.Let us consider the system illustrated by Figure 10, which has been investigated by Awrejcewicz et al.
[41]. Figure 10 shows that a mass , of which the position is denoted as  ∈ R, rests on a conveyor rolling with a constant velocity  ∈ R. The mass is subjected to a nonlinear spring force   () and a rate-dependent friction force   ( ṗ − ).Then, the system is described as the following equation: Here, let us assume that   () and   ( ṗ − ) are defined as follows: where where 0 <  ≪  is a parameter.This approximation was obtained by simply replacing the discontinuity by a linear function of a constant slope in the region || < .Awrejcewicz et al.
[41] has shown that this approximation does reproduce periodic motion appropriately.
Figure 11 shows the simulation results of the proposed approximation (49) and the simple smoothing (50).It shows that the proposed approximation (49) also provides periodic solution, and it is very close to that of the simple smoothing (50).Considering that the result of (50) has been analytically validated through Tikhonov theorem in [41], one can see that the result of the new approximation (49) is also valid.

Conclusion
This paper has introduced a new method to approximate DIs describing nonsmooth mechanical systems involving dry friction and rigid unilateral contact by ODEs.A main difference of the new method from conventional regularization methods is that the resultant ODEs are equivalent to DAIs that are approximations of DIs.As a consequence, the approximated ODEs preserve important features of the original DIs such as static friction and always-repulsive contact force.An algebraic procedure for yielding the ODE approximations has been presented and has been illustrated by using some examples.
Future research should address the theoretical and numerical studies on the influence of the chosen parameters (, , ) on the system behavior.Currently there are no guidelines for the choice of the parameter values; thus they have been chosen through trial and error in the presented examples.In particular, the choice of  and  strongly influences the realized coefficient of restitution.Theoretical or empirical relations between the parameter values and the coefficient of restitution must be sought in the future study.
One limitation of the presented approach is that it is only for "lumped" contacts.In some situations, the contact force may be distributed across a contact area.It is unclear whether the presented approach is applicable or not to such situations.Anisotropic friction force and elastic contact, such as those seen in vehicle tires, would demand further extension of the presented approach.

Figure 1 :
Figure 1: The graphs of relevant functions introduced in Section 2.
Figure7shows the result of the simulation by using(37) with RK4.Here,  1 and  2 are set as high as possible to achieve small penetrations during collisions, and  1 ,  2 , and  are chosen based on some trials and errors.Figures7(a) and 7(b) show bouncing motion in the  direction while Figures 7(c) and 7(d) show a transition from pure translation (slipping in contact) to pure rolling.In Figure7(e), one can find small penetrations produced by the approximation.Moreover, in Figure7(f), one can see impulse oscillations after collisions, which are also consequences of the approximation.Despite these small artifacts, the overall shapes of the graphs in Figures7(a) to 7(d) are close to the expected behaviors of the original DI(33).