# A Differential Algebraic Method to Approximate

Authors: Xiaogang Xiong; Ryo Kikuuwe; Motoji Yamamoto
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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 [5–9] and contact models [10–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–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. Section2 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:sgn, sat, and dio. 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 0 denotes the zero vector of an appropriate dimension.First, let us define the signum functionsgn:ℝn→ℝn and the unit saturation function sat:ℝ+×ℝn→ℝn as follows:
(1)sgn(x)≜{x∥x∥if∥x∥≠0{z∈ℝn∣∥z∥≤1}if∥x∥=0,(2)sat(Z,x)≜{Zx∥x∥if∥x∥>Zxif∥x∥≤Z,
where x∈ℝn, Z∈ℝ+, and ∥·∥ denotes the vector two-norm. If n=1, the sgn(x) and sat(Z,x) 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.The graphs of relevant functions introduced in Section2.
(a)
(b)
(c)
(d)Theorem 1.
Forx,y∈ℝn and Z∈ℝ+, the following relation holds true [27, 28]:
(3)y∈Zsgn(x-y)⟺y=sat(Z,x).Proof.
A proof can be given as follows:(4)y∈Zsgn(x-y)⟺(y=Z(x-y)∥x-y∥∧x≠y)∨(y=x∧∥y∥≤Z)⟺(y=ZxZ+∥x-y∥∧x≠y∧∥y∥=Z)∨(y=x∧∥x∥≤Z)⟺(y=ZxZ+∥x-y∥∧∥x∥=Z+∥x-y∥>Z)∨(y=x∧∥x∥≤Z)⟺(y=Zx∥x∥∧∥x∥>Z)∨(y=x∧∥x∥≤Z)⟺y=sat(Z,x).Next, let us define the “diode” functiondio:ℝ+→ℝ+ as follows:
(5)dio(x)≜{0ifx>0ℝ+ifx=0,
where x∈ℝ+. The following theorem is useful to rewrite DIs involving dio by ODEs.Theorem 2.
Forx∈ℝ and y∈ℝ+, the following relation holds true:
(6)y∈dio(x+y)⟺y=max(0,-x).Proof.
A proof can be given as follows:(7)y∈dio(x+y)⟺(y=0∧x+y>0)∨(y≥0∧x+y=0)⟺(y=0∧x>0)∨(y≥0∧y=-x)⟺y=max(0,-x).The graphs ofdio(x) and max(0,-x) are illustrated as Figures 1(c) and 1(d), respectively.It must be noted that Theorems1 and 2 are special cases of the following relation, which has been used in, for example, [24, AppendixA.3], [29, equation (2)], and [30, equation (4)]:
(8)x-y∈NS(y)⟺y=prox(S,x).
Here, x∈ℝn, y∈S⊂ℝn, S is a closed convex set, NS(y) is the normal cone to the set S at y, and prox(S,x) is the “proximal point” function defined as follows:
(9)prox(S,x)≜argminz∈S∥z-x∥2.
Theorems 1 and 2 can be obtained by using the relation (8) with S={z∈ℝn∣∥z∥≤Z} and S=ℝ+, respectively.


## 3. Previous Approaches

### 3.1. Dry Friction

Let us consider the situation where a rigid massM>0, of which the position is p∈ℝr(r∈{1,2}), is sliding on a fixed surface. Let us assume that it is subject to the dry friction force f∈ℝr and an external force fe∈ℝr. Then, the equation of motion of the mass is described as the following DI:
(10)Mp¨=fe-f,
where
(11)f∈Fsgn(p˙)
and F>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(p˙) is not determined at p˙=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:(12a)a˙=p˙-K∥p˙∥aF,(12b)f=Ka+B(p˙-K∥p˙∥aF)+Dp˙,wherea∈ℝr is the new state variable, K>0 is a sufficiently large constant, and B,D>0 are constants appropriately chosen to suppress the oscillation in p. Dahl friction model [5] is a special case of LuGre friction model with D=B=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–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.


### 3.2. Rigid Unilateral Contact

Let us consider that the one-dimensional system is composed of a rigid massM, of which the position is p∈ℝ, and a fixed rigid wall whose position coincides with the origin. The rigid mass is subject to an external force fe∈ℝ. Then, the equation of motion of the rigid mass is described as the following DI:
(13)Mp¨=fe+f,
where
(14)f∈dio(p).
The integration of (13) and (14) is also difficult due to dio(p), whose value is not determined at p=0.One of the trivial methods to approximately realize the contact forcef in (14) is as follows [24, 34, 35]:
(15)f={-Kp-Bp˙ifp≤00ifp>0,
where K is a sufficiently large positive constant and B is a positive constant large enough to dampen the oscillation in p. This force law can be viewed as a linear viscoelastic contact model with the stiffness K and the viscosity B. 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:
(16)f={max(0,-Kp-Bp˙)ifp≤00ifp>0.
However, both (15) and (16) are discontinuous with respect to p and p˙. 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 p and p˙, but it can also produce unnatural sucking force when p˙ 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.


## 3.1. Dry Friction

Let us consider the situation where a rigid massM>0, of which the position is p∈ℝr(r∈{1,2}), is sliding on a fixed surface. Let us assume that it is subject to the dry friction force f∈ℝr and an external force fe∈ℝr. Then, the equation of motion of the mass is described as the following DI:
(10)Mp¨=fe-f,
where
(11)f∈Fsgn(p˙)
and F>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(p˙) is not determined at p˙=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:(12a)a˙=p˙-K∥p˙∥aF,(12b)f=Ka+B(p˙-K∥p˙∥aF)+Dp˙,wherea∈ℝr is the new state variable, K>0 is a sufficiently large constant, and B,D>0 are constants appropriately chosen to suppress the oscillation in p. Dahl friction model [5] is a special case of LuGre friction model with D=B=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–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.


## 3.2. Rigid Unilateral Contact

Let us consider that the one-dimensional system is composed of a rigid massM, of which the position is p∈ℝ, and a fixed rigid wall whose position coincides with the origin. The rigid mass is subject to an external force fe∈ℝ. Then, the equation of motion of the rigid mass is described as the following DI:
(13)Mp¨=fe+f,
where
(14)f∈dio(p).
The integration of (13) and (14) is also difficult due to dio(p), whose value is not determined at p=0.One of the trivial methods to approximately realize the contact forcef in (14) is as follows [24, 34, 35]:
(15)f={-Kp-Bp˙ifp≤00ifp>0,
where K is a sufficiently large positive constant and B is a positive constant large enough to dampen the oscillation in p. This force law can be viewed as a linear viscoelastic contact model with the stiffness K and the viscosity B. 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:
(16)f={max(0,-Kp-Bp˙)ifp≤00ifp>0.
However, both (15) and (16) are discontinuous with respect to p and p˙. 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 p and p˙, but it can also produce unnatural sucking force when p˙ 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.


## 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:(17a)0∈K(a+βa˙)-Fsgn(p˙-a˙),(17b)f=K(a+βa˙).Here,a∈ℝr is a state variable newly introduced, K>0 is a sufficiently large constant, and β>0 is a constant appropriately chosen to suppress the oscillation in p. A physical interpretation of the approximation ((17a) and (17b)) can be illustrated as Figure 2. A friction force described by Fsgn(p˙-a˙) acts on a massless object whose velocity is p˙-a˙, and a viscoelastic element with the stiffness K and the viscosity Kβ produces the force f in  (17b), which exactly balances the friction force.Figure 2
A physical interpretation of ((17a) and (17b)).In Kikuuwe et al.’s method, (17a) is discretized by Backward-Euler method; for example, a˙ is replaced by (ak-ak-1)/T, where T denotes the timestep size and the subscripts denote time indices, and then it is analytically solved with respect to ak 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:(18a)a˙=(sat(F/K,a+βp˙)-a)β,(18b)f=Ksat(FK,a+βp˙).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:
(19)ddt[p˙pa]=[(fe-Ksat(F/K,a+βp˙))Mp˙(sat(F/K,a+βp˙)-a)β].
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 fe was chosen as
(20)fe={min(0.52,0.8t)Nift<2s0.336+0.144cos(100t)Notherwise,
which is, after t=2 s, oscillatory below the static friction level F=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)).Simulation of the system (10); (a) provided external force fe described as (20); (b) simulation results by RK4 with the timestep size 0.001 s. The parameters are chosen as M=1 kg, F=0.5 N, K=5×103 N/m, and β=2×10-3 s. The initial conditions are p=0 m and p˙=0 m/s.
(a)
(b)It should be mentioned that ((18a) and (18b)) is derived by relaxing the rigid constraint between f and p˙ in (11) by introducing an auxiliary variable a 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:(21a)0∈K(e+βe˙)-dio(p+e)(21b)f=K(e+βe˙),whereK and β are appropriate positive constants and e∈ℝ 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 p+e is connected to the mass through a viscoelastic element with the stiffness K and the viscosity Kβ. Due to the contact, the contact force dio(p+e) acts on the massless object and it balances the force f 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 e˙ cannot be obtained explicitly.Figure 4
A physical interpretation of ((21a) and (21b)).The new approach presented here is to add another termαe˙ to the argument of dio(·), which yields the following DAI:
(22)0∈K(e+βe˙)-dio(p+e+αe˙),(23)f=K(e+βe˙),
where α>0 is another appropriate constant. By using Theorem 2, ((22) and (23)) can be equivalently rewritten as the following ODE:
(24)e˙=max(-eβ,-p+eα),(25)f=Kmax(0,e-β(p+e)α).
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 p, p˙, and e, 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–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 e~≜e+αe˙. By using e~, ((22) and (23)) can be rewritten as follows:(26a)0∈ℒ-1[ℒ[K(e~+βe~˙)](1+αs)-dio(p+e~)],(26b)f=ℒ-1[ℒ[K(e~+βe~˙)]1+αs],whereℒ denotes the Laplace transform. By noting the similarity between ((26a) and (26b)) and ((21a) and (21b)), one can see that force f 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 β=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 state-space model, the system (13) and (14) is approximated by the following ODE:
(27)ddt[p˙pe]=[(fe+Kmax(0,e-β(p+e)/α))Mp˙max(-eβ,-p+eα)].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.Figure 5
Simulation of the system (27) by RK4 with the timestep size 0.001 s. The parameters are chosen as M=1 kg, fe=-9.8 N, K=105 N/m, β=0.01 s, and α=0.01 s (gray dashed), 0.007 s (black dashed), 0.005 s (gray solid), and 0.001 s (black solid). The initial conditions are p=1 m, and p˙=0 m/s.


### 4.3. Dry-Frictional, Rigid Unilateral Contact

The methods in Sections4.1 and 4.2 can be easily combined to describe a rigid unilateral contact involving dry friction. Let us consider a rigid mass M of which the position is p∈ℝ3 and a rigid frictional surface perpendicular to the z-axis and including the origin. Then, the state-space model of the system can be described as the following DI:
(28)ddt[p……p˙]∈[p˙1M[-μdio(pz)sgn(p˙xy)dio(pz)]….……………………………………………],
where p=[pxyT,pz]T, pxy∈ℝ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 (F in (11) ) while that of dio(·) can be done independently. In conclusion, the DI (28) can be approximated by the following ODE:
(29)ddt[p….p˙….ea]=[p˙……………………………………….……………...1M[-fxy(pz,p˙xy,e,a)fz(pz,e)]……………………………………….……………...max(-eβ1,-(pz+e)α)(fxy(pz,p˙xy,e,a)/K2-a)β2],
where
(30)fz(pz,e)≜K1max(0,e-β1(pz+e)α),fxy(pz,p˙xy,e,a)≜K2sat(μfz(pz,e)K2,β2p˙xy+a),
and the parameters α,β1,β2,K1, and K2 are appropriate positive constants. This ODE is obtained by replacing dio(pz) in (28) by fz(pz,e) and then replacing fz(pz,e)sgn(p˙xy) by fxy(pz,p˙xy,e,a).


### 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:(31)x˙∈Φ(x),
where x∈ℝn 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. Let us assume that, in Φ, m different arguments, denoted as ψi(x), i∈{1,…,m}, are used for dio(·) and that l different arguments, denoted as θi(x), i∈{1,…,l}, are used for sgn(·). Here, ψi:ℝn→ℝ+ and θi:ℝn→ℝr, r∈{1,2}, are continuous functions.By applying the methods introduced in Sections4.1 and 4.2, Φ(x) can be approximated by the following procedure.(1)
First, replacedio(ψi(x)) by Kdimax(0,ei-βdi(ψi(x)+ei)/αi), where Kdi, βdi, and αi are appropriate positive constants.(2)
Next, letχi(x) denote the multiplicative factors of sgn(θi(x)), which are nonnegative continuous functions. Then, replace χi(x)sgn(θi) by Ksisat(χi(x)/Ksi,βsiθi(x)+ai), where Ksi and βsi are appropriate positive constants.(3)
Finally, appende˙i=max(-ei/βdi,-(ψi(x)+ei)/αi), i∈{1,…,m}, and a˙i=(sat(χi(x)/Ksi,ai+βsiθi(x))-ai)/βsi, i∈{1,…,l}, to the state-space model.With this procedure, the nonsmooth system (31) is approximated by the following ODE:
(32)ddt[xe1⋮ema1⋮al]=[Φ^(x,e1,…,em,a1,…,al)max(-e1βd1,-ψ1(x)+e1α1)⋮max(-emβdm,-ψm(x)+emαm)(sat(χ1(x)/Ks1,a1+βs1θ1(x))-a1)βs1⋮(sat(χl(x)/Ksl,al+βslθl(x))-al)βsl],
where Φ^(x,e1,…,em,a1,…,al) is the function Φ(x) in which the aforementioned replacements are made.The presented procedure cannot apply if the functionΦ(x) includes a sgn(·) whose multiplicative factor involves discontinuous functions other than dio(·) and if χi(x) are not guaranteed to be nonnegative. The authors, however, are not aware of nonsmooth mechanical systems that must be described by such Φ(x) functions.


## 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:(17a)0∈K(a+βa˙)-Fsgn(p˙-a˙),(17b)f=K(a+βa˙).Here,a∈ℝr is a state variable newly introduced, K>0 is a sufficiently large constant, and β>0 is a constant appropriately chosen to suppress the oscillation in p. A physical interpretation of the approximation ((17a) and (17b)) can be illustrated as Figure 2. A friction force described by Fsgn(p˙-a˙) acts on a massless object whose velocity is p˙-a˙, and a viscoelastic element with the stiffness K and the viscosity Kβ produces the force f in  (17b), which exactly balances the friction force.Figure 2
A physical interpretation of ((17a) and (17b)).In Kikuuwe et al.’s method, (17a) is discretized by Backward-Euler method; for example, a˙ is replaced by (ak-ak-1)/T, where T denotes the timestep size and the subscripts denote time indices, and then it is analytically solved with respect to ak 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:(18a)a˙=(sat(F/K,a+βp˙)-a)β,(18b)f=Ksat(FK,a+βp˙).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:
(19)ddt[p˙pa]=[(fe-Ksat(F/K,a+βp˙))Mp˙(sat(F/K,a+βp˙)-a)β].
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 fe was chosen as
(20)fe={min(0.52,0.8t)Nift<2s0.336+0.144cos(100t)Notherwise,
which is, after t=2 s, oscillatory below the static friction level F=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)).Simulation of the system (10); (a) provided external force fe described as (20); (b) simulation results by RK4 with the timestep size 0.001 s. The parameters are chosen as M=1 kg, F=0.5 N, K=5×103 N/m, and β=2×10-3 s. The initial conditions are p=0 m and p˙=0 m/s.
(a)
(b)It should be mentioned that ((18a) and (18b)) is derived by relaxing the rigid constraint between f and p˙ in (11) by introducing an auxiliary variable a 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:(21a)0∈K(e+βe˙)-dio(p+e)(21b)f=K(e+βe˙),whereK and β are appropriate positive constants and e∈ℝ 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 p+e is connected to the mass through a viscoelastic element with the stiffness K and the viscosity Kβ. Due to the contact, the contact force dio(p+e) acts on the massless object and it balances the force f 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 e˙ cannot be obtained explicitly.Figure 4
A physical interpretation of ((21a) and (21b)).The new approach presented here is to add another termαe˙ to the argument of dio(·), which yields the following DAI:
(22)0∈K(e+βe˙)-dio(p+e+αe˙),(23)f=K(e+βe˙),
where α>0 is another appropriate constant. By using Theorem 2, ((22) and (23)) can be equivalently rewritten as the following ODE:
(24)e˙=max(-eβ,-p+eα),(25)f=Kmax(0,e-β(p+e)α).
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 p, p˙, and e, 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–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 e~≜e+αe˙. By using e~, ((22) and (23)) can be rewritten as follows:(26a)0∈ℒ-1[ℒ[K(e~+βe~˙)](1+αs)-dio(p+e~)],(26b)f=ℒ-1[ℒ[K(e~+βe~˙)]1+αs],whereℒ denotes the Laplace transform. By noting the similarity between ((26a) and (26b)) and ((21a) and (21b)), one can see that force f 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 β=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 state-space model, the system (13) and (14) is approximated by the following ODE:
(27)ddt[p˙pe]=[(fe+Kmax(0,e-β(p+e)/α))Mp˙max(-eβ,-p+eα)].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.Figure 5
Simulation of the system (27) by RK4 with the timestep size 0.001 s. The parameters are chosen as M=1 kg, fe=-9.8 N, K=105 N/m, β=0.01 s, and α=0.01 s (gray dashed), 0.007 s (black dashed), 0.005 s (gray solid), and 0.001 s (black solid). The initial conditions are p=1 m, and p˙=0 m/s.


## 4.3. Dry-Frictional, Rigid Unilateral Contact

The methods in Sections4.1 and 4.2 can be easily combined to describe a rigid unilateral contact involving dry friction. Let us consider a rigid mass M of which the position is p∈ℝ3 and a rigid frictional surface perpendicular to the z-axis and including the origin. Then, the state-space model of the system can be described as the following DI:
(28)ddt[p……p˙]∈[p˙1M[-μdio(pz)sgn(p˙xy)dio(pz)]….……………………………………………],
where p=[pxyT,pz]T, pxy∈ℝ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 (F in (11) ) while that of dio(·) can be done independently. In conclusion, the DI (28) can be approximated by the following ODE:
(29)ddt[p….p˙….ea]=[p˙……………………………………….……………...1M[-fxy(pz,p˙xy,e,a)fz(pz,e)]……………………………………….……………...max(-eβ1,-(pz+e)α)(fxy(pz,p˙xy,e,a)/K2-a)β2],
where
(30)fz(pz,e)≜K1max(0,e-β1(pz+e)α),fxy(pz,p˙xy,e,a)≜K2sat(μfz(pz,e)K2,β2p˙xy+a),
and the parameters α,β1,β2,K1, and K2 are appropriate positive constants. This ODE is obtained by replacing dio(pz) in (28) by fz(pz,e) and then replacing fz(pz,e)sgn(p˙xy) by fxy(pz,p˙xy,e,a).


## 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:(31)x˙∈Φ(x),
where x∈ℝn 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. Let us assume that, in Φ, m different arguments, denoted as ψi(x), i∈{1,…,m}, are used for dio(·) and that l different arguments, denoted as θi(x), i∈{1,…,l}, are used for sgn(·). Here, ψi:ℝn→ℝ+ and θi:ℝn→ℝr, r∈{1,2}, are continuous functions.By applying the methods introduced in Sections4.1 and 4.2, Φ(x) can be approximated by the following procedure.(1)
First, replacedio(ψi(x)) by Kdimax(0,ei-βdi(ψi(x)+ei)/αi), where Kdi, βdi, and αi are appropriate positive constants.(2)
Next, letχi(x) denote the multiplicative factors of sgn(θi(x)), which are nonnegative continuous functions. Then, replace χi(x)sgn(θi) by Ksisat(χi(x)/Ksi,βsiθi(x)+ai), where Ksi and βsi are appropriate positive constants.(3)
Finally, appende˙i=max(-ei/βdi,-(ψi(x)+ei)/αi), i∈{1,…,m}, and a˙i=(sat(χi(x)/Ksi,ai+βsiθi(x))-ai)/βsi, i∈{1,…,l}, to the state-space model.With this procedure, the nonsmooth system (31) is approximated by the following ODE:
(32)ddt[xe1⋮ema1⋮al]=[Φ^(x,e1,…,em,a1,…,al)max(-e1βd1,-ψ1(x)+e1α1)⋮max(-emβdm,-ψm(x)+emαm)(sat(χ1(x)/Ks1,a1+βs1θ1(x))-a1)βs1⋮(sat(χl(x)/Ksl,al+βslθl(x))-al)βsl],
where Φ^(x,e1,…,em,a1,…,al) is the function Φ(x) in which the aforementioned replacements are made.The presented procedure cannot apply if the functionΦ(x) includes a sgn(·) whose multiplicative factor involves discontinuous functions other than dio(·) and if χi(x) are not guaranteed to be nonnegative. The authors, however, are not aware of nonsmooth mechanical systems that must be described by such Φ(x) 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 Figure6. The surface includes the origin and is perpendicular to the z-axis. This example is also introduced in [34] and a similar example is employed in [11]. Let p∈ℝ3 be the position of the gravity center of the object, q be the unit quaternion representing the attitude of the object, and ω∈ℝ3 be the angular velocity of the object. Let R and M 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:
(33)ddt[p˙……....p……….ωq…....….]∈[1M[-μdio(pz-R)sgn(v(p˙,ω))dio(pz-R)]-gp˙.……..………………………………………………………………………………….52MR2(d×[-μdio(pz-R)sgn(v(p˙,ω))dio(pz-R)])…………………………………………………………………………………………Q(ω,q)……….…………………………………………………………………………………],
where
(34)v(p˙,ω)≜[p˙x-ωyR,p˙y+ωxR]T,Q:ℝ3×ℝ4→ℝ4 denotes an appropriate function that transforms ω into the quaternion rate q˙, d≜[0,0,-R]T, and g≜[0,0,9.8]Tm/s2.Figure 6
Example I: a rolling sphere with collision and slip. In the simulation, the parameters are chosen asM=1 kg, R=0.5 m, and μ=0.1 and the initial conditions are p˙=[5.5,0,0]T m/s, p=[0,0,2R]T, and ω=[0,0,0]T rad/s.According to the procedure presented in Section4.4, the DI (33) can be approximated by an ODE in the following procedure. First, one should replace dio(pz-R) by
(35)ψ(pz,e)≜K1max(0,e-β1(e+pz-R)α),
where e∈ℝ and K1, β1, and α are appropriate positive constants. Then μdio(pz-R)sgn(v(p˙,ω)) becomes μψ(pz,e)sgn(v(p˙,ω)). Next, μψ(pz,e)sgn(v(p˙,ω)) should be replaced by
(36)θ(pz,v(p˙,ω),e,a)≜K2sat(μψ(pz,e)K2,a+β2v(p˙,ω)),
where a∈ℝ2 and K2 and β2 are positive constants appropriately chosen. Finally, ODEs defining the behaviors of e and a should be appended to (33). Then, (33) is approximated by the following ODE:
(37)ddt[p˙…...p…....ω……qea]=[1M[-θ(pz,v(p˙,ω),e,a)ψ(pz,e)]-gp˙………………………………………………………………………...52MR2(d×[-θ(pz,v(p˙,ω),e,a)ψ(pz,e)])………………………………………………………………………...Q(ω,q)………………………………………………………………………...max(-eβ1,-e+pz-Rα)(θ(pz,v(p˙,ω),e,a)/K2-a)β2].Figure7 shows the result of the simulation by using (37) with RK4. Here, K1 and K2 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) and 7(b) show bouncing motion in the z direction while Figures 7(c) and 7(d) show a transition from pure translation (slipping in contact) to pure rolling.Figure 7
Simulation results of Example I by using (37) integrated by RK4 with the timestep size 0.001 s. The parameters are chosen as K1=K2=1×105 N/m, β1=β2=4×10-3 s, and α=2.8×10-3 s. Graphs (e) and (f) are enlarged views of graphs (a) and (d), respectively.In Figure7(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 Figure8, which consists of a conveyor moving at a constant velocity u, a spring with the stiffness Ks, two rigid objects M1 and M2, and a rigid vertical wall. The object M1 can move freely in the horizontal direction and is subject to the elastic force from a spring Ks in the vertical direction. It is assumed that the objects do not rotate. The coefficients of friction between the wall and M1, between M1 and M2, and between M2 and the conveyor are μ1, μ2, and μ3, respectively. The state vector is defined as x≜[pT,p˙T]T∈ℝ8 where p=[p1x,p1y,p2x,p2y]T in which [p1x,p1y]T and [p2x,p2y]T denote the positions of M1 and M2, respectively. The object M1 is regarded as being at [0,0]T when it is in contact with the wall and the spring balances the gravity. The object M2 is regarded as being at [0,0]T when it is in contact with the conveyor and the object M1 being at its [0,0]T. Then, the state-space model of the system can be described as the following DI:
(38)ddt[p…....p˙]∈[p˙dio(p1x)-dio(p2x-p1x)M1……………………………………………………….-Ω1(x)-Ω2(x)-Ks1p1yM1-Ω3(x,u)+dio(p2x-p1x)M2Ω2(x)+dio(p2y)M2-g],
where Ω1(x)≜μ1dio(p1x)sgn(p˙1y),Ω2(x)≜μ2dio(p2x-p1x)sgn(p˙1y-p˙2y), and Ω3(x,u)≜μ3dio(p2y)sgn(p˙2x+u).Figure 8
Example II: multiple frictional-unilateral contacts. In the simulation, the parameters were chosen asμ1=μ2=μ3=0.5, M1=0.5 kg, M2=1kg, Ks=100 N/m, and u=1 m/s, and the initial conditions are p=[0,0.25,0.05,0]T m, and p˙=[0,0,0,0]T m/s.According to the procedure presented in Section4.4, the DI (38) is approximated by an ODE in the following procedure. First, one should replace dio(p1x), dio(p2x-p1x), and dio(p2y) by
(39)ψ1(x,e1)≜K1max(0,e1-β1(p1x+e1)α1),(40)ψ2(x,e2)≜K2max(0,e2-β2(p2x-p1x+e2)α2),(41)ψ3(x,e3)≜K3max(0,e3-β3(p2y+e3)α3),
respectively, where Ki, βi, and αi(i∈{1,2,3}) are appropriate positive constants. Then, Ω1(x), Ω2(x), and Ω3(x,u) are found to be replaced by μ1ψ1(x,e1)sgn(p˙1y), μ2ψ2(x,e2)sgn(p˙1y-p˙2y), and μ3ψ3(x,e3)sgn(p˙2x+u), respectively. Next, they should be replaced by θ1(x,e1,a1), θ2(x,e2,a2), and θ3(x,u,e3,a3), respectively, where
(42)θ1(x,e1,a1)≜K4sat(μ1ψ1(x,e1)K4,a1+β4p˙1y),(43)θ2(x,e2,a2)≜K5sat(μ2ψ2(x,e2)K5,a2+β5(p˙1y-p˙2y)),(44)θ3(x,u,e3,a3)≜K6sat(μ3ψ3(x,e3)K6,a3+β6(p˙2x+u)),
and Ki and βi(i∈{4,5,6}) are appropriate constants. Finally, the differential equations defining the behaviors of ei and ai(i∈{1,2,3}) should be appended to (38). Then, (38) is approximated by the following ODE:
(45)ddt[p…...p˙…...e1e2e3a1a2a3]=[p˙ψ1(x,e1)-ψ2(x,e2)M1…………………………………………………………...………………………-Ksp1y-θ1(x,e1,a1)-θ2(x,e2,a2)M1-θ3(x,u,e3,a3)+ψ2(x,e2)M2θ2(x,e2,a2)+ψ3(x,e3)M2-gmax(-e1β1,-(p1x+e1)α1)…………………………………………………………...………………………max(-e2β2,-(p2x-p1x+e2)α2)max(-e3β3,-(p2y+e3)α3)θ1(x,e1,a1)/K4-a1β4θ2(x,e2,a2)/K5-a2β5θ3(x,e3,a3)/K6-a3β6].A numerical simulation was performed by using the ODE (45) with RK4. The results are shown in Figure 9. Here, again, Ki are set as high as possible to achieve small penetrations during collisions, and βi and αi are chosen based on some trials and errors. The time periods indicated by the gray regions are those in which the objects M1 and M2 are in contact to each other. Figures 9(a) and 9(b) show the horizontal bouncing motion of M1 and M2, which eventually converges. Figure 9(c) shows the vertical motion of M1, which exhibits nonsmooth changes in the velocity during the contact with M2, being influenced by the friction force. The vertical position of M1 does not converge to zero because of the static friction forces from the wall and M2. Figure 9(d) shows the vertical motion of M2, which determines the normal force from the conveyor to M2. It properly shows the influence on the normal force from the friction force acting on the side face.Simulation results of Example II by using (45) integrated by RK4 with timestep size 0.001 s. The gray regions indicate the time periods in which the objects M1 and M2 are in contact with each other. The parameters are chosen as Ki=5×105 N/m, βi=2×10-3s(∀i∈{1,…,6}), and αi=1.6×10-3s(∀i∈{1,2,3}).
(a)
(b)
(c)
(d)Also in this simulation, one can observe small penetrations at the time of collisions in Figures9(a) and 9(d). In addition, in Figure 9(b), one can find some impulsive responses in p˙1x, 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 Figure10, which has been investigated by Awrejcewicz et al. [41]. Figure 10 shows that a mass M, of which the position is denoted as p∈ℝ, rests on a conveyor rolling with a constant velocity u∈ℝ. The mass is subjected to a nonlinear spring force Fs(p) and a rate-dependent friction force Fc(p˙-u). Then, the system is described as the following equation:
(46)Mp¨+Fs(p)-Fc(p˙-u)=0.
Here, let us assume that Fs(p) and Fc(p˙-u) are defined as follows:
(47)Fs(p)≜-k1p+k2p3,Fc(p˙-u)≜-FV+|p˙-u|sgn(p˙-u),
where k1, k2, F, and V are positive constants. By using a new variable w≜u-p˙, one can rewrite (46) into the following DI:
(48)ddt[p…...w]∈[u-w1M(k2p3-k1p-FV+|w|sgn(w))………………………………………………………………………].Figure 10
Example III: periodic motion. In the simulation, the parameters were chosen asF=0.2 N, M=1 kg, k1=1 N/m, k2=1N/m3, u=1 m/s, and V=1 m/s, and the initial conditions are p=1.19149 m, and w=0 m/s.According to the procedure presented in Section4.4, (48) is approximated as follows:
(49)ddt[p…...w…...a]=[u-w1M(k2p3-k1p-sat(FV+|w|,Ka+Kβw))…………………………………………………………………………………………...1βsat(F(V+|w|)K,a+βw)-aβ…………………………………………………………………………………………....],
where a is a new state variable and K and β are positive parameters. In contrast, Awrejcewicz et al. [41] used the following equation to approximate (48):
(50)ddt[p…...w]=[u-w1M(k2p3-k1p-Fsat(ϵ,w)(V+max(|w|,ϵ))ϵ)……………………………………………………………………………….],
where 0<ϵ≪V is a parameter. This approximation was obtained by simply replacing the discontinuity by a linear function of a constant slope in the region |w|<ϵ. Awrejcewicz et al. [41] has shown that this approximation does reproduce periodic motion appropriately.Figure11 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.Simulation results of Example III by using (49) and (50) integrated by RK4 with timestep size 0.001 s. The parameters in (49) are chosen as K=1×105 N/m and β=0.5 s. The parameter in (50) is chosen as ϵ=10-5 m/s.
(a)
(b)


## 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 Figure6. The surface includes the origin and is perpendicular to the z-axis. This example is also introduced in [34] and a similar example is employed in [11]. Let p∈ℝ3 be the position of the gravity center of the object, q be the unit quaternion representing the attitude of the object, and ω∈ℝ3 be the angular velocity of the object. Let R and M 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:
(33)ddt[p˙……....p……….ωq…....….]∈[1M[-μdio(pz-R)sgn(v(p˙,ω))dio(pz-R)]-gp˙.……..………………………………………………………………………………….52MR2(d×[-μdio(pz-R)sgn(v(p˙,ω))dio(pz-R)])…………………………………………………………………………………………Q(ω,q)……….…………………………………………………………………………………],
where
(34)v(p˙,ω)≜[p˙x-ωyR,p˙y+ωxR]T,Q:ℝ3×ℝ4→ℝ4 denotes an appropriate function that transforms ω into the quaternion rate q˙, d≜[0,0,-R]T, and g≜[0,0,9.8]Tm/s2.Figure 6
Example I: a rolling sphere with collision and slip. In the simulation, the parameters are chosen asM=1 kg, R=0.5 m, and μ=0.1 and the initial conditions are p˙=[5.5,0,0]T m/s, p=[0,0,2R]T, and ω=[0,0,0]T rad/s.According to the procedure presented in Section4.4, the DI (33) can be approximated by an ODE in the following procedure. First, one should replace dio(pz-R) by
(35)ψ(pz,e)≜K1max(0,e-β1(e+pz-R)α),
where e∈ℝ and K1, β1, and α are appropriate positive constants. Then μdio(pz-R)sgn(v(p˙,ω)) becomes μψ(pz,e)sgn(v(p˙,ω)). Next, μψ(pz,e)sgn(v(p˙,ω)) should be replaced by
(36)θ(pz,v(p˙,ω),e,a)≜K2sat(μψ(pz,e)K2,a+β2v(p˙,ω)),
where a∈ℝ2 and K2 and β2 are positive constants appropriately chosen. Finally, ODEs defining the behaviors of e and a should be appended to (33). Then, (33) is approximated by the following ODE:
(37)ddt[p˙…...p…....ω……qea]=[1M[-θ(pz,v(p˙,ω),e,a)ψ(pz,e)]-gp˙………………………………………………………………………...52MR2(d×[-θ(pz,v(p˙,ω),e,a)ψ(pz,e)])………………………………………………………………………...Q(ω,q)………………………………………………………………………...max(-eβ1,-e+pz-Rα)(θ(pz,v(p˙,ω),e,a)/K2-a)β2].Figure7 shows the result of the simulation by using (37) with RK4. Here, K1 and K2 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) and 7(b) show bouncing motion in the z direction while Figures 7(c) and 7(d) show a transition from pure translation (slipping in contact) to pure rolling.Figure 7
Simulation results of Example I by using (37) integrated by RK4 with the timestep size 0.001 s. The parameters are chosen as K1=K2=1×105 N/m, β1=β2=4×10-3 s, and α=2.8×10-3 s. Graphs (e) and (f) are enlarged views of graphs (a) and (d), respectively.In Figure7(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 Figure8, which consists of a conveyor moving at a constant velocity u, a spring with the stiffness Ks, two rigid objects M1 and M2, and a rigid vertical wall. The object M1 can move freely in the horizontal direction and is subject to the elastic force from a spring Ks in the vertical direction. It is assumed that the objects do not rotate. The coefficients of friction between the wall and M1, between M1 and M2, and between M2 and the conveyor are μ1, μ2, and μ3, respectively. The state vector is defined as x≜[pT,p˙T]T∈ℝ8 where p=[p1x,p1y,p2x,p2y]T in which [p1x,p1y]T and [p2x,p2y]T denote the positions of M1 and M2, respectively. The object M1 is regarded as being at [0,0]T when it is in contact with the wall and the spring balances the gravity. The object M2 is regarded as being at [0,0]T when it is in contact with the conveyor and the object M1 being at its [0,0]T. Then, the state-space model of the system can be described as the following DI:
(38)ddt[p…....p˙]∈[p˙dio(p1x)-dio(p2x-p1x)M1……………………………………………………….-Ω1(x)-Ω2(x)-Ks1p1yM1-Ω3(x,u)+dio(p2x-p1x)M2Ω2(x)+dio(p2y)M2-g],
where Ω1(x)≜μ1dio(p1x)sgn(p˙1y),Ω2(x)≜μ2dio(p2x-p1x)sgn(p˙1y-p˙2y), and Ω3(x,u)≜μ3dio(p2y)sgn(p˙2x+u).Figure 8
Example II: multiple frictional-unilateral contacts. In the simulation, the parameters were chosen asμ1=μ2=μ3=0.5, M1=0.5 kg, M2=1kg, Ks=100 N/m, and u=1 m/s, and the initial conditions are p=[0,0.25,0.05,0]T m, and p˙=[0,0,0,0]T m/s.According to the procedure presented in Section4.4, the DI (38) is approximated by an ODE in the following procedure. First, one should replace dio(p1x), dio(p2x-p1x), and dio(p2y) by
(39)ψ1(x,e1)≜K1max(0,e1-β1(p1x+e1)α1),(40)ψ2(x,e2)≜K2max(0,e2-β2(p2x-p1x+e2)α2),(41)ψ3(x,e3)≜K3max(0,e3-β3(p2y+e3)α3),
respectively, where Ki, βi, and αi(i∈{1,2,3}) are appropriate positive constants. Then, Ω1(x), Ω2(x), and Ω3(x,u) are found to be replaced by μ1ψ1(x,e1)sgn(p˙1y), μ2ψ2(x,e2)sgn(p˙1y-p˙2y), and μ3ψ3(x,e3)sgn(p˙2x+u), respectively. Next, they should be replaced by θ1(x,e1,a1), θ2(x,e2,a2), and θ3(x,u,e3,a3), respectively, where
(42)θ1(x,e1,a1)≜K4sat(μ1ψ1(x,e1)K4,a1+β4p˙1y),(43)θ2(x,e2,a2)≜K5sat(μ2ψ2(x,e2)K5,a2+β5(p˙1y-p˙2y)),(44)θ3(x,u,e3,a3)≜K6sat(μ3ψ3(x,e3)K6,a3+β6(p˙2x+u)),
and Ki and βi(i∈{4,5,6}) are appropriate constants. Finally, the differential equations defining the behaviors of ei and ai(i∈{1,2,3}) should be appended to (38). Then, (38) is approximated by the following ODE:
(45)ddt[p…...p˙…...e1e2e3a1a2a3]=[p˙ψ1(x,e1)-ψ2(x,e2)M1…………………………………………………………...………………………-Ksp1y-θ1(x,e1,a1)-θ2(x,e2,a2)M1-θ3(x,u,e3,a3)+ψ2(x,e2)M2θ2(x,e2,a2)+ψ3(x,e3)M2-gmax(-e1β1,-(p1x+e1)α1)…………………………………………………………...………………………max(-e2β2,-(p2x-p1x+e2)α2)max(-e3β3,-(p2y+e3)α3)θ1(x,e1,a1)/K4-a1β4θ2(x,e2,a2)/K5-a2β5θ3(x,e3,a3)/K6-a3β6].A numerical simulation was performed by using the ODE (45) with RK4. The results are shown in Figure 9. Here, again, Ki are set as high as possible to achieve small penetrations during collisions, and βi and αi are chosen based on some trials and errors. The time periods indicated by the gray regions are those in which the objects M1 and M2 are in contact to each other. Figures 9(a) and 9(b) show the horizontal bouncing motion of M1 and M2, which eventually converges. Figure 9(c) shows the vertical motion of M1, which exhibits nonsmooth changes in the velocity during the contact with M2, being influenced by the friction force. The vertical position of M1 does not converge to zero because of the static friction forces from the wall and M2. Figure 9(d) shows the vertical motion of M2, which determines the normal force from the conveyor to M2. It properly shows the influence on the normal force from the friction force acting on the side face.Simulation results of Example II by using (45) integrated by RK4 with timestep size 0.001 s. The gray regions indicate the time periods in which the objects M1 and M2 are in contact with each other. The parameters are chosen as Ki=5×105 N/m, βi=2×10-3s(∀i∈{1,…,6}), and αi=1.6×10-3s(∀i∈{1,2,3}).
(a)
(b)
(c)
(d)Also in this simulation, one can observe small penetrations at the time of collisions in Figures9(a) and 9(d). In addition, in Figure 9(b), one can find some impulsive responses in p˙1x, 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 Figure10, which has been investigated by Awrejcewicz et al. [41]. Figure 10 shows that a mass M, of which the position is denoted as p∈ℝ, rests on a conveyor rolling with a constant velocity u∈ℝ. The mass is subjected to a nonlinear spring force Fs(p) and a rate-dependent friction force Fc(p˙-u). Then, the system is described as the following equation:
(46)Mp¨+Fs(p)-Fc(p˙-u)=0.
Here, let us assume that Fs(p) and Fc(p˙-u) are defined as follows:
(47)Fs(p)≜-k1p+k2p3,Fc(p˙-u)≜-FV+|p˙-u|sgn(p˙-u),
where k1, k2, F, and V are positive constants. By using a new variable w≜u-p˙, one can rewrite (46) into the following DI:
(48)ddt[p…...w]∈[u-w1M(k2p3-k1p-FV+|w|sgn(w))………………………………………………………………………].Figure 10
Example III: periodic motion. In the simulation, the parameters were chosen asF=0.2 N, M=1 kg, k1=1 N/m, k2=1N/m3, u=1 m/s, and V=1 m/s, and the initial conditions are p=1.19149 m, and w=0 m/s.According to the procedure presented in Section4.4, (48) is approximated as follows:
(49)ddt[p…...w…...a]=[u-w1M(k2p3-k1p-sat(FV+|w|,Ka+Kβw))…………………………………………………………………………………………...1βsat(F(V+|w|)K,a+βw)-aβ…………………………………………………………………………………………....],
where a is a new state variable and K and β are positive parameters. In contrast, Awrejcewicz et al. [41] used the following equation to approximate (48):
(50)ddt[p…...w]=[u-w1M(k2p3-k1p-Fsat(ϵ,w)(V+max(|w|,ϵ))ϵ)……………………………………………………………………………….],
where 0<ϵ≪V is a parameter. This approximation was obtained by simply replacing the discontinuity by a linear function of a constant slope in the region |w|<ϵ. Awrejcewicz et al. [41] has shown that this approximation does reproduce periodic motion appropriately.Figure11 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.Simulation results of Example III by using (49) and (50) integrated by RK4 with timestep size 0.001 s. The parameters in (49) are chosen as K=1×105 N/m and β=0.5 s. The parameter in (50) is chosen as ϵ=10-5 m/s.
(a)
(b)


## 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 (K, α, β) 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.



---

*Source: 320276-2013-05-08.xml*