Abstract

Nonsmooth mechanical systems, which are mechanical systems involving dry friction and rigid unilateral contact, are usually described as differential inclusions (DIs), that is, differential equations involving discontinuities. Those DIs may be approximated by ordinary differential equations (ODEs) by simply smoothing the discontinuities. Such approximations, however, can produce unrealistic behaviors because the discontinuous natures of the original DIs are lost. This paper presents a new algebraic procedure to approximate DIs describing nonsmooth mechanical systems by ODEs with preserving the discontinuities. The procedure is based on the fact that the DIs can be approximated by differential algebraic inclusions (DAIs), and thus they can be equivalently rewritten as ODEs. The procedure is illustrated by some examples of nonsmooth mechanical systems with simulation results obtained by the fourth-order Runge-Kutta method.

1. 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 penalty-based 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 [59] and contact models [1013] 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 [1422] or analytically ([23, III.A], [24, 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.

2. Mathematical Preliminaries

For the discussion throughout this paper, this section introduces three functions: , , and . Some theorems regarding those functions are also presented. In the rest of this paper, denotes the set of all real numbers and denotes the set of all nonnegative real numbers. The symbol denotes the zero vector of an appropriate dimension.

First, let us define the signum function and the unit saturation function as follows: where , , and denotes the vector two-norm. If , the and can be depicted as Figures 1(a) and 1(b), respectively. The following theorem is useful to rewrite the DIs involving as ODEs involving sat.

Theorem 1. For and , the following relation holds true [27, 28]:

Proof. A proof can be given as follows:

Next, let us define the “diode” function as follows: where . The following theorem is useful to rewrite DIs involving dio by ODEs.

Theorem 2. For and , the following relation holds true:

Proof. A proof can be given as follows:

The graphs of and 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, A.3], [29, equation ], and [30, equation ]: Here, , , is a closed convex set, is the normal cone to the set at , and is the “proximal point” function defined as follows: Theorems 1 and 2 can be obtained by using the relation (8) with and , respectively.

3. Previous Approaches

3.1. Dry Friction

Let us consider the situation where a rigid mass , of which the position is , is sliding on a fixed surface. Let us assume that it is subject to the dry friction force and an external force . Then, the equation of motion of the mass is described as the following DI: where and 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 is not determined at , according to the definition (1) of .

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 is the new state variable, is a sufficiently large constant, and are constants appropriately chosen to suppress the oscillation in . Dahl friction model [5] is a special case of LuGre friction model with . 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, 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, 1719]. One exception is Kikuuwe et al.'s approach [23, 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.

3.2. Rigid Unilateral Contact

Let us consider that the one-dimensional system is composed of a rigid mass , of which the position is , and a fixed rigid wall whose position coincides with the origin. The rigid mass is subject to an external force . 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 , whose value is not determined at .

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, 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.

4. 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.

4.1. 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, is a state variable newly introduced, is a sufficiently large constant, and 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 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.

In Kikuuwe et al.’s method, (17a) is discretized by Backward-Euler method; for example, is replaced by , 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  s, oscillatory below the static friction level  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)).

4.2. 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 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 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 , which yields the following DAI: where 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 [1113], 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 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 low-pass filtered viscoelastic force although it does not exist in the real world. When , ((26a) and (26b)) is equivalent to ((21a) and (21b)) with , 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 state-space 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.

4.3. Dry-Frictional, Rigid Unilateral Contact

The methods in Sections 4.1 and 4.2 can be easily combined to describe a rigid unilateral contact involving dry friction. Let us consider a rigid mass of which the position is and a rigid frictional surface perpendicular to the -axis and including the origin. Then, the state-space model of the system can be described as the following DI: where , , and is the friction coefficient between the mass and the surface.

It must be noticed that (28) includes a multiplication of and . To approximate this, one must replace first and then replace because the replacement of involves its multiplicative factor ( in (11) ) while that of can be done independently. In conclusion, the DI (28) can be approximated by the following ODE: where and the parameters , and are appropriate positive constants. This ODE is obtained by replacing in (28) by and then replacing by .

4.4. 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 is the state vector of the system. Here, is a function that contains and in several places and may also contain single valued functions. Let us assume that, in , different arguments, denoted as , , are used for and that different arguments, denoted as , , are used for . Here, and , , are continuous functions.

By applying the methods introduced in Sections 4.1 and 4.2, can be approximated by the following procedure.(1)First, replace by , where , , and are appropriate positive constants.(2)Next, let denote the multiplicative factors of , which are nonnegative continuous functions. Then, replace by , where and are appropriate positive constants.(3)Finally, append , , and , , to the state-space model.

With this procedure, the nonsmooth system (31) is approximated by the following ODE: where is the function in which the aforementioned replacements are made.

The presented procedure cannot apply if the function includes a whose multiplicative factor involves discontinuous functions other than 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.

5. Examples

5.1. 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 be the position of the gravity center of the object, be the unit quaternion representing the attitude of the object, and 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: where denotes an appropriate function that transforms into the quaternion rate , , and .

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 by where and , , and are appropriate positive constants. Then becomes . Next, should be replaced by where and and 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, and are set as high as possible to achieve small penetrations during collisions, and , , and are chosen based on some trials and errors. Figures 7(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 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).

5.2. 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 and , and a rigid vertical wall. The object 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 , between and , and between and the conveyor are , , and , respectively. The state vector is defined as where in which and denote the positions of and , respectively. The object is regarded as being at when it is in contact with the wall and the spring balances the gravity. The object is regarded as being at when it is in contact with the conveyor and the object being at its . Then, the state-space model of the system can be described as the following DI: where , and .

According to the procedure presented in Section 4.4, the DI (38) is approximated by an ODE in the following procedure. First, one should replace , , and by respectively, where , , and are appropriate positive constants. Then, , , and are found to be replaced by , , and , respectively. Next, they should be replaced by , , and , respectively, where and and are appropriate constants. Finally, the differential equations defining the behaviors of and 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 and are in contact to each other. Figures 9(a) and 9(b) show the horizontal bouncing motion of and , which eventually converges. Figure 9(c) shows the vertical motion of , which exhibits nonsmooth changes in the velocity during the contact with , being influenced by the friction force. The vertical position of does not converge to zero because of the static friction forces from the wall and . Figure 9(d) shows the vertical motion of , which determines the normal force from the conveyor to . 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 , 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).

5.3. 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 , rests on a conveyor rolling with a constant velocity . 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 , , , and are positive constants. By using a new variable , one can rewrite (46) into the following DI:

According to the procedure presented in Section 4.4, (48) is approximated as follows: where is a new state variable and and are positive parameters. In contrast, Awrejcewicz et al. [41] used the following equation to approximate (48): where 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.

6. 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.

Acknowledgment

This work was supported by Grant-in-Aid for Scientific Research B (no. 24360098) from Japan Society for the Promotion of Science (JSPS).