Dynamic simulation of task constrained of a rigid-flexible manipulator

Abstract

A rigid-flexible manipulator may be assigned tasks in a moving environment where the winds or vibrations affect the position and/or orientation of surface of operation. Consequently, losses of the contact and perhaps degradation of the performance may occur as references are changed. When the environment is moving, knowledge of the angle α between the contact surface and the horizontal is required at every instant. In this paper, different profiles for the time varying angle α are proposed to investigate the effect of this change into the contact force and the joint torques of a rigid-flexible manipulator. The coefficients of the equation of the proposed rotating surface are changing with time to determine the new X and Y coordinates of the moving surface as the surface rotates.

1. Introduction

In the last two decades, robotic constrained motion applications, where the end-effector of the robot is always in contact with its environment, have received a considerable attention. Typical examples of such applications include grinding, deburring, cutting, polishing, drilling, fastening, etc. To execute these tasks successfully it is necessary to control both motion of the robot and the contact force between the end-effector and the environment. Despite the voluminous research done in the last decade, the study of dynamics and control of the constrained motion of the flexible manipulators remains open for further investigations.

Modeling, simulation, and control of manipulators in constrained motion attract many researchers for the last two decades. In 1981, the first attempt in dealing with this problem for rigid manipulators was proposed by applying hybrid position/forcce control scheme (Raibert, M. & Craig, J., 1981). For rigid-flexible manipulator, the work by (Su, et al. 1990) represented the earliest study in force control of constrained maneuver. In 1994, a dynamic modeling and force/motion control of constrained flexible robot was presented (Ulsoy, G. & Hu, F., 1994). In their model, the flexible link has an additional degree-of-freedom to move axially into the first rigid link. The dynamic hybrid position/force control of a rigid-flexible (Matsuno, F. & Yamamoto, K., 1994) and two-link flexible manipulator (Matsuno, et al, 1994) where the elastic deformations were approximated by B-spline functions were Investigated. The inverse force and motion control of constrained three- axis elastic robot based on non-linear inversion and stabilization was presented (Yim, W. & Singh, H., 1995). Using a two-time scale force/position controller, a general systemic model of flexible robot interacting with a rigid environment was addressed (Rocco, P. & Book, W., 1996). The recent work by (Shi, et al, 1999)], A mathematical model of a constrained rigid-flexible manipulator based on Hamilton's principle was derived. A multi variable controller was proposed for the simultaneous motion and force control was also investigated. A parallel force and position control scheme of flexible manipulators using perturbation theory was proposed (Siciliano, B. & Villani, L., 2000). In 2000, an algorithm employing the computed torque method in the free space and hybrid motion/force controller for rigid manipulators was presented (Bennon, et al., 2000). Experimental results on non-moving and moving environment were also illustrated. While in 2001, a position and force control scheme for a class of flexible joint robots during constrained motion tasks with model uncertainties based on singular perturbation theory was developed (Hu, Y. & Vukovich, G. 2001).

Dynamic modeling and simulation of constrained motion of rigid-flexible manipulator in contact with a compliant surface was addressed (Ata, A. & and Ghazy, S., 2001). While in 2003, they continued their work towards finding an optimal motion trajectory for the constrained motion based on the minimum energy consumption (Ata, et al., 2003). In this study, an extension to the results in the previous work is presented in such away to further study the interaction of the joint motion profile and the angular motion profile of the constraint surface. The constrained surface is assumed to be rotating, therefore, different angular profiles for the surface are proposed to investigate the effect of this change into the constrained force and the joint torques.

2. Dynamic Modeling

Consider the two-link planar manipulator shown in Figure 1. The first link is rigid and the second link is assumed to be flexible. The end-effector of the robot is in contact with a moving environment during the task. The tip payload consists of mass m₃ with moment of inertia I₃ about its own axis of rotation. The flexible link is assumed to be thin and slender so it can be modeled as an Euler-Bernoulli's beam of length l₂, uniform mass density ρ₂, the cross section I₂. Applying the Virtual Link Coordinate System (VLCS), consider θ₂(t) be the angle of a line pointing from the second joint to the tip mass as described in reference (Benati, M. & Morro, A., 1994). The transversal bending deflection w(x,t) at a point xε[0, l₂] along the second link is described with respect to the virtual link. Longitudinal deformations are neglected. No damping is assumed and the manipulator moves in the horizontal plane so, the gravity is not considered.

Fig. 1:

A two link rigid-flexible manipulator with tip mass

The equations of motion can be derived using the extended Hamilton's principle (Ata, et al., 2003). The constrained motion of the rigid-flexible can be described as:

$$\tau =\mathbf{M}\mathbf{(}\theta \mathbf{)}\ddot{\theta }\mathbf{+}\mathbf{V}\mathbf{(}\theta \mathbf{,}\stackrel{\mathbf{̇}}{\theta }\mathbf{)}\mathbf{+}{\tau }_{\mathbf{r}}\mathbf{+}{\tau }_{\mathbf{f}}$$

(1)

where:

$$\mathbf{\tau }=\left[\begin{array}{c}{\mathbf{\tau }}_{\mathbf{1}}\\ {\mathbf{\tau }}_{\mathbf{2}}\end{array}\right]$$

in which τ₁ and τ₂ are the rigid and flexible torques respectively.

M(θ) is a 2 by 2 inertia matrix of the manipulator and its elements can be described by:

$$\begin{array}{l}{\mathbf{M}}_{\mathbf{11}}=({\mathbf{I}}_{\mathbf{1}}+{\mathbf{I}}_{\mathbf{h1}}+{\mathbf{m}}_{\mathbf{2}}{\mathbf{l}}_{\mathbf{1}}^{\mathbf{2}}+{\mathbf{m}}_{\mathbf{3}}{\mathbf{l}}_{\mathbf{1}}^{\mathbf{2}})\\ {\mathbf{M}}_{\mathbf{12}}=({\mathbf{m}}_{\mathbf{2}}{\mathbf{l}}_{\mathbf{1}}{\bar{\mathbf{x}}}_{\mathbf{2}}+{\mathbf{m}}_{\mathbf{3}}{\mathbf{l}}_{\mathbf{1}}{\mathbf{l}}_{\mathbf{2}})\mathbf{cos}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})\\ -\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\int }}\mathbf{\rho }{\mathbf{l}}_{\mathbf{1}}\mathbf{w}\mathbf{sin}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})\mathbf{\text{dx}}\\ {\mathbf{M}}_{\mathbf{21}}={\mathbf{M}}_{\mathbf{12}}\\ {\mathbf{M}}_{\mathbf{22}}=({\mathbf{I}}_{\mathbf{2}}+{\mathbf{I}}_{\mathbf{h}\mathbf{2}}+{\mathbf{m}}_{\mathbf{3}}{\mathbf{l}}_{\mathbf{2}}^{\mathbf{2}})+\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\int }}\mathbf{\rho }{\mathbf{w}}^{\mathbf{2}}\mathbf{\text{dx}}\end{array}$$

$\mathbf{V}\mathbf{(}\theta \mathbf{,}\stackrel{\mathbf{̇}}{\theta }\mathbf{)}$ is the Coriolis and centripetal terms and can be given by the column vector:

$$\begin{array}{l}\mathbf{V}\left(\mathbf{\theta },\dot{\mathbf{\theta }}\right)=\left(\begin{array}{c}{\mathbf{v}}_{\mathbf{11}}\\ {\mathbf{v}}_{\mathbf{21}}\end{array}\right)\\ {\mathbf{v}}_{\mathbf{11}}=-({\mathbf{m}}_{\mathbf{2}}{\mathbf{l}}_{\mathbf{1}}{\bar{\mathbf{x}}}_{\mathbf{2}}+{\mathbf{m}}_{\mathbf{3}}{\mathbf{l}}_{\mathbf{1}}{\mathbf{l}}_{\mathbf{2}}){\dot{\mathbf{\theta }}}_{\mathbf{2}}^{\mathbf{2}}\mathbf{sin}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})\\ +\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\int }}\begin{array}{c}{\mathbf{\rho }}_{\mathbf{2}}{\mathbf{l}}_{\mathbf{1}}[(\ddot{\mathbf{w}}-\mathbf{w}{\dot{\mathbf{\theta }}}_{\mathbf{2}}^{\mathbf{2}})\mathbf{cos}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})\\ -(\mathbf{2}\dot{\mathbf{w}}{\dot{\mathbf{\theta }}}_{\mathbf{2}})\mathbf{sin}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})]\mathbf{\text{dx}}\end{array}\\ {\mathbf{v}}_{\mathbf{21}}=({\mathbf{m}}_{\mathbf{2}}{\mathbf{l}}_{\mathbf{1}}{\bar{\mathbf{x}}}_{\mathbf{2}}+{\mathbf{m}}_{\mathbf{3}}{\mathbf{l}}_{\mathbf{1}}{\mathbf{l}}_{\mathbf{2}}){\dot{\mathbf{\theta }}}_{\mathbf{1}}^{\mathbf{2}}\mathbf{sin}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})\\ +\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\int }}{\mathbf{\rho }}_{\mathbf{2}}\{-\mathbf{2}\mathbf{w}\dot{\mathbf{w}}{\dot{\mathbf{\theta }}}_{\mathbf{2}}+\mathbf{x}\ddot{\mathbf{w}}+{\mathbf{l}}_{\mathbf{1}}\mathbf{w}{\dot{\mathbf{\theta }}}_{\mathbf{1}}^{\mathbf{2}}\mathbf{cos}({\mathbf{\theta }}_{\mathbf{2}}-{\mathbf{\theta }}_{\mathbf{1}})\}\mathbf{\text{dx}}\end{array}$$

τ_f is the torque due to friction. Frictional force is a highly nonlinear phenomenon that is difficult to model accurately. In the present study, we neglect the effect of this friction force.

τ_r is the interaction torque due to the contact with the environment and given as:

$${\mathbf{\tau }}_{\mathbf{r}}={\mathbf{J}}^{\mathbf{T}}\mathbf{F}={\mathbf{J}}^{\mathbf{T}}\left[\begin{array}{c}{\mathbf{F}}_{\mathbf{x}}\\ {\mathbf{F}}_{\mathbf{y}}\end{array}\right]$$

where J^T is the transpose Jacobian of the system.

$${\mathbf{J}}^{\mathbf{T}}=\left(\begin{array}{l}-{\mathbf{I}}_{\mathbf{1}}\mathbf{sin}{\mathbf{\theta }}_{\mathbf{1}}-{\mathbf{l}}_{\mathbf{2}}\mathbf{sin}{\mathbf{\theta }}_{\mathbf{2}}\hspace{1em}{\mathbf{l}}_{\mathbf{1}}\mathbf{cos}{\mathbf{\theta }}_{\mathbf{1}}+{\mathbf{l}}_{\mathbf{2}}\mathbf{cos}{\mathbf{\theta }}_{\mathbf{2}}\hfill \\ -{\mathbf{l}}_{\mathbf{2}}\mathbf{sin}{\mathbf{\theta }}_{\mathbf{2}}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\mathbf{l}}_{\mathbf{2}}\mathbf{cos}{\mathbf{\theta }}_{\mathbf{2}}\hfill \end{array}\right)$$

F_x and F_y represent the applied forces by the end-effector on the surface in contact including the inertia forces of the end-effector and is modeled with a spring (Su, et al., 1990):

$$\begin{array}{rl}& {\mathbf{F}}_{\mathbf{x}}\mathbf{=}{\mathbf{K}}_{\mathbf{s}}\mathbf{(}{\mathbf{r}}_{\mathbf{x}}\mathbf{-}{\mathbf{l}}_{\mathbf{1}}\mathbf{cos}{\mathit{\theta }}_{\mathbf{1}}\mathbf{-}{\mathbf{l}}_{\mathbf{2}}\mathbf{cos}{\mathit{\theta }}_{\mathbf{2}}\mathbf{)}\mathbf{+}\\ & {\mathbf{m}}_{\mathbf{3}}\mathbf{(}{\mathbf{l}}_{\mathbf{1}}{\ddot{\mathit{\theta }}}_{\mathbf{1}}\mathbf{sin}{\mathit{\theta }}_{\mathbf{1}}\mathbf{+}{\mathbf{l}}_{\mathbf{1}}{\stackrel{\mathbf{̇}}{\mathit{\theta }}}_{\mathbf{1}}^{\mathbf{2}}\mathbf{cos}{\mathit{\theta }}_{\mathbf{1}}\mathbf{+}{\mathbf{l}}_{\mathbf{2}}{\ddot{\mathit{\theta }}}_{\mathbf{2}}\mathbf{sin}{\mathit{\theta }}_{\mathbf{2}}\\ & \mathbf{+}{\mathbf{l}}_{\mathbf{2}}{\stackrel{\mathbf{̇}}{\mathit{\theta }}}_{\mathbf{2}}^{\mathbf{2}}\mathbf{cos}{\mathit{\theta }}_{\mathbf{2}}\mathbf{)}\end{array}$$

(2)

$$\begin{array}{rl}& {\mathbf{F}}_{\mathbf{y}}\mathbf{=}{\mathbf{K}}_{\mathbf{s}}\mathbf{(}{\mathbf{r}}_{\mathbf{y}}\mathbf{-}{\mathbf{l}}_{\mathbf{1}}\mathbf{sin}{\mathit{\theta }}_{\mathbf{1}}\mathbf{-}{\mathbf{l}}_{\mathbf{2}}\mathbf{sin}{\mathit{\theta }}_{\mathbf{2}}\mathbf{)}\\ & \mathbf{+}{\mathbf{m}}_{\mathbf{3}}\mathbf{(}{\mathbf{l}}_{\mathbf{1}}{\stackrel{\mathbf{̇}}{\mathit{\theta }}}_{\mathbf{1}}^{\mathbf{2}}\mathbf{sin}{\mathit{\theta }}_{\mathbf{1}}\mathbf{-}{\mathbf{l}}_{\mathbf{1}}{\ddot{\mathit{\theta }}}_{\mathbf{1}}\mathbf{cos}{\mathit{\theta }}_{\mathbf{1}}\\ & \mathbf{+}{\mathbf{l}}_{\mathbf{2}}{\stackrel{\mathbf{̇}}{\mathit{\theta }}}_{\mathbf{2}}^{\mathbf{2}}\mathbf{sin}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathbf{l}}_{\mathbf{2}}\ddot{{\mathit{\theta }}_{\mathbf{2}}}\mathbf{cos}{\mathit{\theta }}_{\mathbf{2}}\mathbf{)}\end{array}$$

(3)

Where K_s is the spring stiffness (N/m), r_x, and r_y are the coordinates of the contact point.

For the second link, the equation due to the flexibility effect is given by:

$$\begin{array}{rl}& {\mathit{\rho }}_{\mathbf{2}}\{\mathbf{w}{\stackrel{\mathbf{̇}}{\mathit{\theta }}}_{\mathbf{2}}^{\mathbf{2}}\mathbf{-}\ddot{\mathbf{w}}\mathbf{-}\mathbf{x}{\ddot{\mathit{\theta }}}_{\mathbf{2}}\mathbf{-}{\mathbf{l}}_{\mathbf{1}}\mathbf{[}\stackrel{\mathbf{̇}}{{\mathit{\theta }}_{\mathbf{1}}^{\mathbf{2}}}\mathbf{sin}\mathbf{(}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathit{\theta }}_{\mathbf{1}}\mathbf{)}\\ & \mathbf{+}{\ddot{\mathit{\theta }}}_{\mathbf{1}}\mathbf{cos}\mathbf{(}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathit{\theta }}_{\mathbf{1}}\mathbf{)}\mathbf{]}\mathbf{\}}\mathbf{-}{\mathbf{k}}_{\mathbf{2}}\mathbf{w}\mathbf{″}\mathbf{″}\mathbf{=}\mathbf{0}\end{array}$$

(4)

subject to four boundary conditions.

$$\mathbf{w}\mathbf{(}\mathbf{0}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}\mathbf{0}$$

(5a)

$$\mathbf{w}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}\mathbf{0}$$

(5b)

$${\mathbf{w}}^{\mathbf{\prime }\mathbf{\prime }}\mathbf{(}\mathbf{0}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}\mathbf{-}\frac{{\mathit{\tau }}_{\mathbf{21}}\mathbf{(}\mathbf{t}\mathbf{)}}{{\mathbf{k}}_{\mathbf{2}}}$$

(5c)

$${\mathbf{w}}^{\mathbf{\prime }\mathbf{\prime }}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}{\mathit{\omega }}^{\mathbf{2}}{\mathbf{v}}^{\mathbf{\prime }}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{)}{\mathbf{I}}_{\mathbf{3}}\mathbf{(}\mathbf{t}\mathbf{)}$$

(5d)

The fourth boundary condition is due to the existence of the inertia of the tip mass around its own axis of rotation (Clough, R. & Penzien, J., 1993) in which, ω is the natural frequency of the flexible link, v′(l₂) is an assumed function of spatial coordinate, and I₃(t) is the inertia of tip mass which is a time function because we consider tip mass varies linearly with time.

3. Vibration of the flexible link with time-dependent boundary conditions

For an angular displacement θ(t) and an elastic deflection w(x,t), the total displacement y(x,t) of a point along the flexible link from the hub can be considered as:

$$\mathbf{y}\mathbf{(}\mathbf{x}\mathbf{,}\mathbf{t}\mathbf{)}=\mathbf{w}\mathbf{(}\mathbf{x}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{+}\mathbf{x}{\mathit{\theta }}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}$$

(6)

By ignoring the first term of (4) since its effect is only obvious at very high speed, and substituting from equation (6) into equations (4 and 5) one can get:

$$\begin{array}{rl}& {\mathit{\rho }}_{\mathbf{2}}\mathbf{\{}\mathbf{-}\ddot{\mathbf{y}}\mathbf{-}{\mathbf{l}}_{\mathbf{1}}\mathbf{[}{\dot{\mathit{\theta }}}_{\mathbf{1}}^{\mathbf{2}}\mathbf{sin}\mathbf{(}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathit{\theta }}_{\mathbf{1}}\mathbf{)}\\ & \mathbf{+}{\ddot{\mathit{\theta }}}_{\mathbf{1}}\mathbf{cos}\mathbf{(}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathit{\theta }}_{\mathbf{1}}\mathbf{)}\mathbf{]}\mathbf{\}}\mathbf{-}{\mathbf{E}}_{\mathbf{2}}{\mathbf{I}}_{\mathbf{2}}{\mathbf{y}}^{\mathbf{\prime }\mathbf{\prime }\mathbf{\prime }\mathbf{\prime }}\mathbf{=}\mathbf{0}\end{array}$$

(7)

$$\mathbf{y}\mathbf{(}\mathbf{0}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}\mathbf{0}$$

(8a)

$$\mathbf{y}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}{\mathbf{l}}_{\mathbf{2}}{\mathit{\theta }}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}$$

(8b)

$${\mathbf{y}}^{\mathbf{\prime }\mathbf{\prime }}\mathbf{(}\mathbf{0}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}\mathbf{-}\frac{{\mathit{\tau }}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}}{{\mathbf{E}}_{\mathbf{2}}{\mathbf{I}}_{\mathbf{2}}}$$

(8c)

$${\mathbf{y}}^{\mathbf{\prime }\mathbf{\prime }}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{,}\mathbf{t}\mathbf{)}\mathbf{=}{\mathit{\omega }}^{\mathbf{2}}{\mathbf{v}}^{\mathbf{\prime }}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{)}{\mathbf{I}}_{\mathbf{3}}\mathbf{(}\mathbf{t}\mathbf{)}$$

(8d)

Our objective is to find the rigid and flexible hub torques to move the end-effector through a prescribed trajectory by solving the inverse dynamics problem. Solving equations (1 and 4) for the joints torque subject to the boundary conditions (5a-d) to obtain the rigid and flexible hub torques is a very difficult task. This is simply because one has to calculate the elastic deflection of the arm to obtain the flexible hub torque. Unfortunately, the required flexible torque is also included in the time-dependent boundary condition (8c). An alternative approach to the computation of the link deformation is to use approximations for the flexible torque (Asada, et al. 1990). The sequence of calculations to get τ₁ and τ₂, can be summarized as follows:

Assume the joint motion profiles θ₁(t) and θ₂ (t).

Consider τ₂ as a rigid torque without any elastic effect, substitute into equation (8c).

Solve equation (7) subject to the boundary conditions (8) to get y(x,t) and then substitute into (6) to get w(x,t), $\dot{\mathbf{w}}(\mathbf{x}\mathbf{,}\mathbf{t}\mathbf{)}$ and $\ddot{\mathbf{w}}(\mathbf{x},\mathbf{t})$

Calculate F_x and F_y from equations (2 & 3) and substitute into (1) to get τ₁ and τ₂.

For the rigid-flexible manipulator under consideration, spatial discretization techniques (e. g., assumed modes, finite element or Galerkin's method) are typically employed to obtain a finite dimensional system of ordinary differential equation suitable for simulation (Hu and Ulsoy, 1994). The assumed modes method, based on modal expansion, consists in finding a particular solution (admissible function) of the Euler-Bernoulli equation to satisfy only the geometric boundary conditions. On the other hand, the finite element method consists in finding a local solution of the equation over a finite element of the flexible link. In this research we are applying the assumed modes method as a technique for the analytical solution.

The solution for y(x,t) can be obtained using the assumed modes method [Meirovitch, L., 1967 & Low, K., 1989] in the form:

$$\begin{array}{rl}& {\mathbf{y}}_{\mathbf{n}}(\mathbf{x},\mathbf{t})=\sum _{\mathbf{0}}^{\mathbf{n}}[{\mathbf{v}}_{\mathbf{n}}(\mathbf{x}){\zeta }_{\mathbf{n}}(\mathbf{t})]+\mathbf{g}(\mathbf{x})\mathbf{e}(\mathbf{t})\\ & \hspace{0.17em}+\mathbf{h}(\mathbf{x})\mathbf{f}(\mathbf{t})+\mathbf{q}(\mathbf{x})\mathbf{p}(\mathbf{t})\end{array}$$

(9)

$$\mathbf{e}\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{=}{\mathit{\theta }}_{\mathbf{2}}\mathbf{(}\mathbf{t}\mathbf{)}$$

(10a)

$$\mathbf{f}\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{=}\mathbf{-}\frac{{\mathit{\tau }}_{\mathbf{2}}}{{\mathbf{E}}_{\mathbf{2}}{\mathbf{I}}_{\mathbf{2}}}$$

(10b)

$$\mathbf{p}\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{=}{\mathit{\omega }}^{\mathbf{2}}{\mathbf{v}}^{\mathbf{\prime }}\mathbf{(}{\mathbf{l}}_{\mathbf{2}}\mathbf{)}{\mathbf{l}}_{\mathbf{3}}(\mathrm{t})$$

(10c)

The subscript n indicates that the continuous system has been approximated by an n-degree-of-freedom system. The three functions g (x), h (x) and q (x) are of the spatial coordinate alone to satisfy the homogeneous boundary conditions for v_n (x), and ζ_n (t) is the time function. Taking the lowest possible order of polynomials the functions g (x), h (x) and q (x) can be written as:

$$\mathbf{g}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{=}\mathbf{x}$$

(11a)

$$\mathbf{h}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{=}\mathbf{-}\frac{\mathbf{1}}{\mathbf{3}}{\mathbf{l}}_{\mathbf{2}}\mathbf{x}\mathbf{+}\frac{\mathbf{1}}{\mathbf{2}}{\mathbf{x}}^{\mathbf{2}}\mathbf{-}{\frac{\mathbf{1}}{\mathbf{61}}}_{\mathbf{2}}{\mathbf{x}}^{\mathbf{3}}$$

(11b)

$$\mathbf{q}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{=}\mathbf{-}\frac{\mathbf{1}}{\mathbf{6}}{\mathbf{l}}_{\mathbf{2}}\mathbf{x}\mathbf{+}{\frac{\mathbf{1}}{\mathbf{61}}}_{\mathbf{2}}{\mathbf{x}}^{\mathbf{3}}$$

(11c)

The corresponding eigenfunction v_n (x) can be obtained as [14]:

$${\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{=}\sqrt{\frac{\mathbf{2}}{\mathit{\rho }{\mathbf{l}}_{\mathbf{2}}}}\mathbf{sin}\frac{\mathbf{n}\mathit{\pi }}{{\mathbf{l}}_{\mathbf{2}}}\mathbf{x}\mathbf{,}\mathbf{n}\mathbf{=}\mathbf{1}\mathbf{,}\mathbf{2}\mathbf{,}\mathbf{3}\mathbf{\dots }$$

(12)

The last two nonlinear terms inside the parentheses in equation (6) can be regarded as distributed excitation force with unit density. This effect can be compensated in the time function (Meirovitch, L., 1967) as:

$${\mathit{\zeta }}_{\mathbf{n}}\mathbf{=}\frac{\mathbf{1}}{{\mathit{\omega }}_{\mathbf{n}}}\underset{\mathbf{0}}{\overset{\mathbf{t}}{\mathbf{\int }}}{\mathbf{N}}_{\mathbf{n}}\mathbf{(}\mathit{\tau }\mathbf{)}\mathbf{sin}\mathbf{(}\mathbf{t}\mathbf{-}\mathit{\tau }\mathbf{)}\mathbf{d}\mathit{\tau }$$

(13)

where

$$\begin{array}{rl}& {\mathbf{N}}_{\mathbf{n}}\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{=}{\mathbf{N}}_{\mathbf{n}\mathbf{l}}\mathbf{(}\mathbf{t}\mathbf{)}\mathbf{-}{\mathit{\rho }}_{\mathbf{2}}{\mathbf{l}}_{\mathbf{1}}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}\\ & \mathbf{[}{\ddot{\mathit{\theta }}}_{\mathbf{1}}\mathbf{cos}\mathbf{(}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathit{\theta }}_{\mathbf{1}}\mathbf{)}\mathbf{+}{\stackrel{\mathbf{̇}}{\mathit{\theta }}}_{\mathbf{1}}^{\mathbf{2}}\mathbf{sin}\mathbf{(}{\mathit{\theta }}_{\mathbf{2}}\mathbf{-}{\mathit{\theta }}_{\mathbf{1}}\mathbf{)}\mathbf{]}\mathbf{d}\mathbf{x}\\ & {\mathbf{G}}_{\mathbf{n}}^{\mathbf{\ast }}\mathbf{=}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}{\mathbf{g}}^{\mathbf{\prime }\mathbf{\prime }\mathbf{\prime }\mathbf{\prime }}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{d}\mathbf{x}{\mathbf{G}}_{\mathbf{n}}\mathbf{=}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}{\mathit{\rho }}_{\mathbf{2}}\mathbf{g}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{d}\mathbf{x}\\ & {\mathbf{H}}_{\mathbf{n}}^{\mathbf{\ast }}\mathbf{=}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}{\mathbf{h}}^{\mathbf{\prime }\mathbf{\prime }\mathbf{\prime }}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{d}\mathbf{x}{\mathbf{H}}_{\mathbf{n}}\mathbf{=}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}{\mathit{\rho }}_{\mathbf{2}}\mathbf{h}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{d}\mathbf{x}\\ & {\mathbf{Q}}_{\mathbf{n}}^{\mathbf{\ast }}\mathbf{=}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}{\mathbf{q}}^{\mathbf{\prime }\mathbf{\prime }\mathbf{\prime }}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{d}\mathbf{x}{\mathbf{Q}}_{\mathbf{n}}\mathbf{=}\underset{\mathbf{0}}{\overset{{\mathbf{l}}_{\mathbf{2}}}{\mathbf{\int }}}{\mathbf{v}}_{\mathbf{n}}\mathbf{(}\mathbf{x}\mathbf{)}{\mathit{\rho }}_{\mathbf{2}}\mathbf{q}\mathbf{(}\mathbf{x}\mathbf{)}\mathbf{d}\mathbf{x}\end{array}$$

(14)

The convolution integral (14) can be evaluated using Duhamel integral method (Clough, R. & Penzien, J., 1993). The simulation algorithm thus developed is implemented in Matlab and the results will be presented in the next section.

4. Simulation results and discussion

For the solution of the inverse dynamics problem, three joint motion profiles, namely sine, polynomial, and Gaussian velocity profiles are assumed. In order to investigate the interaction of the two motion profiles on the contact force and the joint torques, four different profiles for the rotation of the constrained surface have been applied to each joint motion profile. These four profiles have in common starting and ending values but they differ in their rates of increase. Accordingly, the angle α of the constrained surface will start from 2π/3 and reach 5π/6 in 5 seconds which is the same time duration as for joint motion. The proposed contact surface is a parabolic shape: y = ax² + bx + c. The coefficients of the proposed surface are changing with time to determine the new X and Y coordinates of the contact point of the moving surface as the surface rotates. Figures 2 and 3 represent the position of the stationary and rotating surface respectively in X-Y coordinates.

Fig. 2

Stationary surface

Fig. 3

Rotating surface

The system parameters used in this study are given in below:

l₁=0.5 m. l₂=0.75 m,

m₁=0.5 Kg, m₂= l₂*ρ₂, m₃=0.15+(0.15)t/T

I₁=0.0834 Kg.m², I₂=(l₂³*ρ₂/3), I₃₌m_3*K² Kg.m²

I_h1=10*I₁, I_h2=10*I₂

E₂I₂=2.4507 Nm²

Where T is the time duration of motion and K is the radius of gyration for the tip mass.

The joint torques and the contact force have been simulated first for the stationary surface in Figs 4, 6, 8 to show that they are greatly different from those resulted in Figs 5(a-d), 7(a-d) and 9(a-d) for the rotating surface with different angular velocity profiles.

Fig. 4

Sinusoidal profile joint motion on stationary surface

Fig. 5a

Sine-Sine

Fig. 5b

Gaussian-Sine

Fig. 5c

Poly-Sine

Fig. 5d

Parabolic-Sine

Fig. 6

Polynomial profile joint motion on stationary surface

Fig. 7a

Sine-Polynomial

Fig. 7b

Gaussian-Poly.

Fig. 7c

Poly-Poly.

Fig. 7d

Parabolic-Poly.

Fig. 8

Gaussian profile joint motion on stationary surface

Fig. 9a

Sine-Gaussian

Fig. 9b

Gauss.-Gauss.

Fig. 9c

Poly-Gauss.

Fig. 9d

Parabolic-Gauss

Simulation plots are presented in groups according to the joint motion profiles. In the Figs 5(a-d), 7(a-d) and 9(a-d), the first word in the Figure's title indicates the angular motion profile of the constrained surface and the second stands for the joint motion profile.

It can be observsed from Figs 5, 7 and 9 that the contact force and the joint torques are certainly affected by the changing of the angular motion profiles for the constrained surface. In all three profiles of joint motion, the contact force tends to increase. While for sinusoidal profile of the constrained surface (Figs. 5a, 7a, 9a), the contact force increases considerably creating uncertainty about the stability the system at the end of the time interval. The polynomial, parabolic and Gaussian profiles of the constrained surface give better distributions for the contact force and joints torque.

On the other hand, the contact force distributions are influenced by the interaction of the angular velocity profile of the environment and the joint motion profile. This can be observed especially in Figs. 9(a-d) where the distributions of the constrained forces are not as smooth as in Figs. 5(a-d) and 7(a-d) due to the Gaussian joint motion profile. As for the joint torques, all four different angular velocity profiles of the constrained surface produce almost same trend with a little bit different peak values in each particular joint motion profile.

In the case of stationary surface, maximum joint torques and contact force result from Gaussian joint motion profile (Fig. 8), followed by sine profile (Fig. 6) and polynomial profile (Fig. 4).

5. Conclusion

In this paper, the effects of different angular velocity profiles for the constrained surface on the constrained force and the joint torques have been investigated. Simulation results show that system performance can also be affected due to the interaction between the joint motion and the angular motion of the constrained surface.

Therefore, for the better dynamic system performance, it is necessary to choose the motion profiles with judgment for both the surface rotation and the joint motion since there are two motions involved. However, knowledge of the time varying angle α is essential at every instant of surface rotation to obtain correct force measurements and thus joint torques.