Sage Journals: Discover world-class research

Abstract

In this study, a biomechanical model mimicking the human hand-arm system under heavy disk excitation is developed to define the stability threshold between the preload vibration and the stress of the human hand-arm system. The fully electromechanical system, consisting of a DC motor, a transmission system, and three discrete masses representing the upper arm, forearm, and hand lifting a heavy disk, is connected by shock absorbers and various springs. The challenge of mimicking the angular activity of the elbow joint in torsion and flexion was also included to assess its contribution to system stability and to study the absorption of mechanical energy in the human hand-arm system. Finally, the movements of the model were described by a differential matrix equation, and its analysis facilitates the understanding of the threshold of the chaotic behaviors observed in an oscillating arms system. Specifically, it explores how the motion of a DC motor can be regulated using a single controller parameter within the context of a multi-degree-of-freedom human hand-arm system. The study uses numerical integrations to analyze system behaviors, with an emphasis on the use of frequency spectra, Poincaré maps, and bifurcation diagrams for visualization and interpretation. The results indicate that the system exhibits chaotic behavior when subjected to certain conditions, particularly when the mass load exceeds 35 kg. Imposing motion on the mass block further intensifies the chaos, leading to manifestations such as strong oscillations and frequency-synchronization effects. These characteristics, exhibited by the idealized model, which imitates the human body through a mechanical model and computation of the motions of the body, play a vital role in various fields of ergonomics and sports biomechanics. Understanding the dynamic behaviors of such systems, especially in response to varying conditions, has significant implications for engineering applications.

Keywords

DC motor articulated arm shoulder bifurcation analysis Poincaré maps chaotic dynamical systems

Introduction

The hand and arm constitute a very complex, continuous, and heterogeneous system comprising skin, muscles, and bones, adapted to various daily activities. To establish robust standards for tolerable vibration stress and exposure duration, it is imperative to improve an understanding of hand-arm system vibrations, particularly when interacting with external tools. This requires a comprehensive stability analysis and examination of the vibration responses of various elements of the hand-arm system. Robotic arms, designed to resemble human arms in structure and operation, are driven by a combination of mechanical and electrical components. They aim to facilitate smooth cooperation between humans and robots by mimicking their movements and functions. Research has focused on creating robotic arms that resemble human arms to accurately mimic human arm movements and gestures. The effectiveness of human-robot arms hinges on their ability to work harmoniously with humans, which underscores their importance and success in various applications. Dynamic modeling is used to elucidate the relationship between force and motion within a system. Such modeling is generally based on two fundamental principles, as explained in the publication.¹ The first principle, the Euler–Lagrange formulation, evaluates the system as a whole considering the overall kinetic and potential energy. On the contrary, the Newton-Euler formulation dissects the system into individual components. The choice of method depends on various factors such as the number and positioning of joints in the manipulator. For example, in Al-Qahtani et al.,² the dynamic modeling of a robotic manipulator with four degrees-of-freedom (d.o.f) is detailed using the Euler-Lagrange formulation, with the base oriented horizontally. Additionally, Yuguang and Fan³ used a combination of the Lagrange equation and the Morison formula to derive a dynamic model for a three d.o.f manipulator suitable for hydrostatic environments. In particular, the modeling process includes accounting for additional forces acting on the joints when the manipulator operates underwater. Human joints, which consist of various tissues such as cartilage, ligaments, and synovial fluid, interact to create internal loads transmitted to adjacent structures. Joints facilitate angular movement, contributing to the d.o.f problem, where humans can choose motor task solutions without compromising efficacy.⁴ Despite their complexity, each human joint can be approximated by a mechanical model with similar functional capabilities. The selection of the appropriate joint model is crucial for valid and reliable results. Models of mechanical joints, both ideal and non-ideal, have been developed, with skeletal models typically using ideal joints (e.g. ball, universal, hinge). However, complex joints, such as the shoulder or elbow, require sophisticated modeling in the form of “closed-loop” systems, considering coupled rotations.^5–7 Human-robot arms have applications in various industries, including manufacturing, assistive technology, and healthcare. Lovasz et al.⁸ emphasize the importance of precise control and haptic feedback in human-robot arm applications, while Aljinovic et al.⁹ illustrate how integrating human-robot arms into assembly lines can improve productivity and reduce worker fatigue. Challenges persist to ensure effective and safe physical interactions, with advances in artificial intelligence and machine learning contributing to the development of context-aware and adaptive control algorithms for human-robot arms. Effective control strategies are essential for precise and safe interactions between human-robot arms. Liao and Wang¹⁰ have explored impedance control techniques, while Buerkle et al.¹¹ have explored the application of human intention recognition algorithms to anticipate and assist with arm movements. Nonlinear methods are increasingly being applied to accurately assess the stability and capture the intricate dynamics of these interactions. Stability analysis is crucial in human-robot arm-robot interactions, as it ensures that the system maintains balance and avoids unwanted behaviors. Non-linear analysis methods are used to evaluate the joint stability of the human-robot system in various interaction scenarios, providing a deeper understanding of the complex interactions between human operators and robotic arms. Real-time stability analysis is essential to ensure the safe operation of human-robot arms. Bascetta and Ferretti¹² employ Lyapunov-based analysis and predictive control algorithms to assess stability in real-time, significantly contributing to operational safety and security. One challenge of stability analysis is considering disturbances and uncertainties, and researchers like Faramin and Ataei¹³ are investigating robust stability analysis methods that consider external disturbances and parameter uncertainties.

Future directions involve the integration of machine learning algorithms to improve stability predictions. The literature on the application of the bifurcation technique to the analysis of human-robot arm instability provides information on the complex behaviors and possible instabilities that can arise in these systems. Understanding the circumstances under which human-robot arms can display unstable behaviors depends on the study of instability. The bifurcation technique examines how the behavior of the system varies with dissimilarities in key parameters, providing a unique perspective. Gierlak and Warmiński¹⁴ used bifurcation analysis to investigate how human operators and robot arms interact under various load scenarios. This method helps to find bifurcation points that suggest a possible instability. Madsen et al.¹⁵ and other studies focus on creating dynamic models that take into account nonlinearity and coupling between human and robotic arms. Bifurcation techniques are then applied to these models to explore the parameter space and discover regions of instability. Park and Choi¹⁶ examine how different human operator behaviors and control inputs affect joint system stability. Bifurcation diagrams that provide information about the transition from stable to unstable regions are revealed by bifurcation analysis. Krishan and Verma¹⁷ use bifurcation techniques to assess possible instability and ensure safe operation to address security concerns. These analyses facilitate the design of control strategies that avoid system instability and improve operational security. Challenges in bifurcation-based instability analysis include dealing with high-dimensional parameter spaces and uncertainties. Tartaruga et al.¹⁸ explore techniques to handle these challenges, including sensitivity analysis and uncertainty quantification. Ongoing directions involve combining bifurcation analysis with Poincaré for more comprehensive instability predictions and aim to refine Poincaré techniques, enhancing our ability to predict and mitigate instabilities in human-robot arm interactions.

In the fields of sports biomechanics and ergonomics, successful emulation of the human body has been achieved using mechanical models to compute body motions, joint movements, and muscle forces. Lee¹⁹ introduced a biomechanical model of the human upper body incorporating bones and muscles. The model aims to model and control joint bones and muscles to simulate physics-based soft tissue deformations, while also incorporating detailed 3D anatomical geometries of hard and soft tissues. Multibody dynamics is a biomechanical technique widely used in sports, clinical settings, rehabilitation, and ergonomics to improve human performance,^20,21 improve the prosthetic design,²² and detect limb asymmetries after injuries,²⁰ and characterize the loads associated with repetitive movements.²³ It is reasonable to pursue simulation of the vibration response of the hand-arm system using similar methodologies. Therefore, the objective of this article is to construct a mathematical model designed to generate Poincaré and bifurcation signals associated with linear and nonlinear vibrational phenomena in human arm mechanisms. The paper makes a significant contribution by taking a comprehensive approach to the modeling and manipulation of dynamic biomechanical arm systems with elastic universal joints to model the elbow, with emphasis on their sensitivity to external stresses. The approach involves the use of a numerical tool to produce these signals, and particular emphasis will be placed on refining the initial parameters to ensure that their representation aligns with signals obtained experimentally in the literature. The subsequent phase of the study involves the use of these signals in purely numerical simulations of the mathematical model. These simulations aim to predict nonlinear vibrational phenomena, including contact or stability problems among coupled arm components.

The following sections of the paper are structured as follows: Section 2 describes the dynamic model of the electrical servomotor device responsible for rotating machine parts coupled to the two arms. Following this, Section 3 comprehensively describes the entire dynamic system and formulates the mathematical equations that govern the behavior of double arms with flexible joints. In Section 4, numerical simulation results are presented for a specific set of parameters to validate the efficacy of the proposed arms model. Finally, Section 5 provides a concise conclusion.

Description of the dynamic model

A descriptive biomechanical human hand-arm model is presented here, where masses represent the hand and the forearm, and the DC motor with the gear system represents the support, which mimics the shoulder of the body. An integrative approach is adopted for the development and validation of coupled arms under variable loads. The article elucidates the dynamics of human arms featuring five d.o.f, with the dynamic model of the manipulator’s arm derived through the application of Lagrangian formalism. This section details model selection and summarizes its essential features.

Biomechanical model description

The human hand-arm system presented in Figure 1(a) is recognized as one of the fundamental kinematic structures, consisting of a planar arm with two rotating joints that move within the X-Y plane. The schematic diagram shown in Figure 1(b) illustrates the kinematic structure of the arm system, presenting its fundamental elements, their spatial arrangement within the overall system, and its key parameters. More precisely, the diagram delimits two rigid arms characterized by lengths denoted L₁ and L₂, each having a mass represented, respectively, by M₁ and M₂. Additionally, the distance between the joint rotation axis and the center of gravity of each arm is designated by a₁ and a₂. The moments of inertia of the arm bodies (J₁ and J₂) are with respect to an axis that passes through the center of gravity and parallel to the horizontal axis X. The deviation of the arm from its initial position (where the joint angle is not equal to 0°) is described by the angles $φ_{1}$ , $φ_{2}$ , and $φ_{3}$ representing the human hand-arm system d.o.f. Neglecting their lateral deflection and friction effects, only the torsional vibrations of the system are studied. Initially, the two shafts and hinge joints are assumed to be in phase, which means that they are initially aligned in a vertical plane.

Figure 1.

(a) Schematic drawing of the hand-arm system, 1. Contracted bicep, 2. Ulna, 3. Radius, 4. Elbow joint, 5. Humerus, 6. Shoulder and (b) corresponding simplified mechanical model.

The elbow (coupling joint) is a hinge joint that allows forearm rotation. The movement produced by the joint describes the compression and extension of the rotating arms, $φ_{1}$ and $φ_{2}$ , meaning the increase or decrease in the angle between the humerus and the forearm around the flexion-extension axis of the elbow through the torsional deflection $ϕ_{1}$ and $ϕ_{2}$ in a concise range of angles from [−90°, 90°]. The wrist, which connects the hand to the forearm and supports the load to be lifted, is modeled in this work as a one-d.o.f $φ_{3}$ spherical joint to simplify the mathematical analyses related to kinematics. It enables a rotation that covers a range of angles from [−12°, 12°] effectively making the elbow a universal joint. This defined biomechanical model has a significant limitation related to external force, which considerably restricts the natural rotation of the wrist.

Electromechanical control system components (DC motor)

The proposed model involves a DC motor connected to a gearbox system that provides torque and speed to move the elbow arm. In this study, the DC motor operates by electromagnetic induction to produce rotational motion. In the closed-loop control system, the motor’s output (such as speed or position) is continuously monitored to provide feedback on motor performance, allowing adjustments to the input voltage or current. Closed-loop control systems offer greater accuracy, reliability, and stability, especially in dynamic or variable operating conditions. Understanding and optimizing control parameters is crucial to achieving desired performance characteristics in a biomechanical arm system. The DC motor and gear system support imitate the shoulder attached to the body. The shoulder (motor) is considered to generate a rotational movement that can be associated with Euler angle notation $θ_{m}$ , typically covering a range of angles from [−90°, 90°] as shown in Figure 2. The rotation of the motor drives the gear train, which moves the hand-arm system. Figure 2(a) shows an equivalent circuit diagram representing the armature as a circuit with resistance (R_a), inductance (L_a), and voltage source (E_b).

Figure 2.

Schematic diagram of a DC machine: (a) motor schematic and (b) transmission gear system.

The armature speed of a DC motor is regulated by the armature voltage e_a = K_te_v, with Kirchhoff’s voltage law allowing the calculation of the total torque. The relationship between armature current $i_{a} (t)$ applied voltage $e_{a} (t)$ , and back EMF, $v_{b} (t)$ is established as:

(R_{a} + L_{a} s) I_{a} (t) + V_{b} (s) = E_{a} (s)

(1)

The expression for the motor torque (T_m) is formulated in Laplace form, taking into account the proportional relationship with the armature current (I_a) and the impact of inductance (L_a), as outlined by Ogata²⁴.

T_{m} (s) = K_{t} (1 - L_{a} / R_{a} s) I_{a} (s)

(2)

This analysis focuses on the use of gears in robotics applications to reduce the load’s angular speed by a factor of n, analyzing output parameters such as internal damping C_meq and gear inertia J_meq, respectively. For this gear-train system displayed in Figure 2(b), the equation for torque equilibrium with respect to $θ_{m}$ can be expressed as:

(J_{meq} s^{2} + c_{meq} s) θ_{m} (s) = K_{t} (1 - L_{a} / R_{a} s) I_{a} (s)

(3)

The desired transfer function can be obtained by simplifying equations (1) and (2) and removing $I_{a} (s)$ and writing a loop equation around the Laplace transforms armature circuit.

\begin{matrix} \frac{θ_{L} (s)}{E_{a} (s)} = n K_{t} (\frac{s}{J_{meq} R_{a}} - \frac{1}{J_{meq} L_{a}}) / \\ s (s^{2} + (\frac{R_{a}}{L_{a}} + \frac{c_{meq}}{J_{meq}} + \frac{K_{t} K_{b}}{J_{meq} R_{a}}) s + \frac{c_{meq} R_{a}}{J_{meq} L_{a}} - \frac{K_{t} K_{b}}{J_{meq} L_{a}}) \end{matrix}

(4)

Where $J_{meq} = J_{m} + n_{1}^{2} J_{G} + n_{1}^{2} n_{2}^{2} J_{A}$ ; $c_{meq} = c_{m} + n_{1}^{2} c_{G} + n_{1}^{2} n_{2}^{2} c_{A}$ ; and $n = n_{1} n_{2} = (N_{1} N_{3} / N_{2} N_{4})$ are respectively the equivalent moment of inertia; the equivalent viscous-friction coefficient of the gear train relative to the shaft of the motor and the gear train ratio product. The transfer function was derived from a linear system and did not involve the process of linearization around operating points; therefore, the closed-loop transfer function is computed using unity feedback, resulting in the following transfer function result.

\frac{θ_{L} (s)}{E_{a} (s)} = \frac{18000 s - 18}{s (s^{2} + (17 + K_{b}) s + 1000 - 2.964 K_{b})}

(5)

In this work, the utility of the transfer function concept emerged from the identification of fluctuations associated with oscillations in the power system that could potentially be attenuated by a flexible mechanical transmission system. Force characteristics in system analysis are crucial to understanding inherent properties such as inertia, damping, and stiffness. Accurate predictive models help reproduce system dynamics, such as the driving force of an arm, as demonstrated by Bernstein²⁵ and expressed by:

F_{b} (t) = \frac{K_{t} N_{g} V_{b} (t)}{R_{a} r_{mg}} - \frac{K_{t} N_{g} K_{b} \overset{\cdot}{X} (t))}{R_{a} r_{mg}^{2}}

(6)

Here, N_g and r_mg are, respectively, the ratio of the planetary gearbox and the radius of the motor pinion.

In this study, a representative set of values of the DC motor parameters is chosen from Table 1.

Table 1.

The values of the physical parameters of the DC motor.

DC motor-gearbox	Value and unit
Motor moment of inertia (J_meq)	0.000005 kg/m²
Motor viscous friction constant (C_meq)	0.000346 kg/m²
Motor torque constant (K_t)	0.00767 N.m/Amp
Electric resistance (R_a)	2.6 Ω
Electric inductance (L_a)	0.5 H
Gear ratio (n)	0.116
Electromotive back EMF (K_b)	Adjustable(V-s/rad)

Implement a feedback mechanism defined in equation (5), it is possible to continuously monitor the state of the system and adjust the control input accordingly by varying K_b. A first simulation allows us to obtain the unit step response of the transfer function given by equation (5), in the unstable cycle, which further represents the angular displacement of the output shaft and is presented in Figure 3. A family of transfer function curves with different values of adjustable K_b is shown in Figure 3(a). The Nyquist diagram associated with K_b shows the convergence of the response graphs of Figure 3(a) and Nyquist Figure 3(b) toward the point of instability. The presence of closed-loop poles in the left half-plane, as shown in Figure 3(b), confirms the stability of the system. This indicates that K_b = 1.92 is the critical value for the appearance of limit cycles. Similarly, on the Nyquist diagram associated with K_b = 80.98, the graph of K_b passes through the point (0, 0), signifying the imminent disappearance of the associated limit cycle. Figure 3 illustrates how the speed of the DC motor can be controlled using adjusted controllers (K_b). The K_b controller offers better performance in both transient and steady-state responses. From the results of the motor speed simulation displayed in Figure 3(c) and (d), the parameter K_b value may result in greater sensitivity to load changes and less stable speed control. The simulation in Figure 3(c) diverges from the original equilibrium point, in agreement with the predictions of the initial investigation. Adjusting the back EMF. The constant K_b can enhance stability, as shown in the unit-step response in Figure 3(a) and Nyquist plots in Figure 3(b), where a symmetric and stable ellipsis indicates linear and stable behavior.

Observation 1: The evolution of the Nyquist plot demonstrated that the system was stable for a defined disturbance, since the polar diagram does not go around the axis in one direction and the poles are in the left half-plane (negative real part). The K_b parameter significantly influences the behavior of a DC motor, affecting the voltage and speed relationships. After some trial and error, K_b = 1982 provides a satisfactory response. Adjusting the K_b parameter is crucial for optimal system performance, with higher values enhancing stability and speed regulation, while lower values may decrease back EMF, causing load sensitivity and less stable speed control.

Figure 3.

The transfer function. It is seen that the third set of parameters K_b controls the system with greater precision: (a) time (s), (b) real axis, (c) time (s), and (d) time (s).

Establishment of the joint disturbance function

For the system considered, although they have limitations compared to more specialized components, the shock absorbers and springs are integrated into biomechanical models to simulate specific aspects of human muscle-tendon-ligament structures. To validate the model, the behavior of the joint is compared to that of a cardan shaft. By performing the inverse Laplace transform, the time-domain representation of the DC motor response, allowing the analysis and interpretation of its behavior over time, can be obtained as $ϕ_{1}$ . The analysis consists of giving a displacement to the input arm; the resulting deformation is then measured according to different angles of rotation of the joint, as shown in Figure 4.

Figure 4.

Simplified modeling of the elbow joint in flexion and extension: (a) arm and bone anatomy, (b) arm in flexion, and (c) arm in extension.

The model incorporates both a torsion spring and a damper, which effectively reproduce the arm flexion in Figure 4(a) and the extension in Figure 4(b). Flexion decreases the angle, while extension increases it, as shown in Figure 4(a) and (b). To simulate different postures of the arm, it is assumed that the joint angle $δ$ is significant, around 30° or more, to meaningfully facilitate the linearization of this model. Given that, the displacement of inertia $J_{b 2}$ is defined by:

ϕ_{2} (t) = ϕ_{1} (t) + ξ ϕ_{1} (t) = ϕ_{1} (t) (1 + ξ (t))

(7)

This study uses $ξ$ as a disturbance parameter, modeled using Cardan joint kinematics, to analyze the relationship between the upper and lower arm Tchomeni Kouejou and Alugongo²⁶.

\tan ϕ_{2} (t) = \cos δ (t) \tan ϕ_{1} (t)

(8)

The equation can be solved by referencing equation (8) and applying the tangent trigonometric function after appropriate manipulation.

\tan (ϕ_{1} (t) ξ (t)) = \frac{(1 + \cos δ (t)) \tan ϕ_{1} (t)}{(\cos δ (t) \tan^{2} ϕ_{1} (t) + 1)}

(9)

The model of a flexible joint arm manipulator, given in equation (10), defines the perturbation function for a significant variation angle. This is determined by rearranging equation (9) and solving for $ξ$ as the subject to obtain Tchomeni Kouejou and Alugongo²⁶:

\begin{matrix} ξ (ϕ_{1}) = \\ {\begin{matrix} 0 if ϕ_{1} = mod π \\ \frac{1}{ϕ_{1} (t)} \tan^{- 1} (\frac{(1 + \cos δ) \tan ϕ_{1} (t)}{(\cos δ \tan^{2} ϕ_{1} (t) + 1)}) if \frac{ϕ_{1}}{π} + \frac{1}{2} \notin Z, ϕ_{1} \neq 0 \end{matrix} \end{matrix}

(10)

The characteristic dynamic behaviors of the disturbance function $ξ$ can finally be analyzed by varying the angle of the joint. Figure 5 shows the fluctuation response of the simulated disturbance joint function where the initial common angle is 20°.

Figure 5.

(a) Arm joint disturbance for four values of $δ$ in flexion and extension positions and (b) bifurcation diagram of the system (10) with changes in the K_b parameter with a fixed value $δ$ .

The angles were established for individual oscillations in the interval [20°; 180°]. Figure 5(a) illustrates the deformations measured according to different angles of rotation of the joint and the bifurcation obtained numerically. The joint oscillates, causing periodic changes and misalignment angle variations. High values of the parameters increase the instability of arm movement, with angles above 60° showing significant instability. The observation aligns with Figure 5(b), which shows a map of the bifurcation diagram of the perturbation function. Applying the first iterative method with a step length of 0.01, equation (9) shows that around $ϕ_{1} (K_{b}) = 0.8$ , the motion becomes increasingly sensitive to initial conditions and eventually stabilizes as K_b approaches 4. Thus, the region near $ϕ_{1} (K_{b}) = 0.95$ , where the system transforms into a weakly nonlinear oscillator is identified as the transition point from chaos to stabilization. In general, arm joint disorders tend to increase with higher values of the misalignment angle, whether the joint is in flexion or extension. This disturbance can compromise stability, depending on the severity of the misalignment. The rotation angle $ϕ_{1}$ is therefore large and very fluctuating.

Chaotic behaviour of the joint disturbance function

To examine the impact of disturbance parameters on system motion, a bifurcation diagram can be generated using MATLAB. Figure 6 shows that the joint disturbance function given by equation (10) where $0 \leq δ \leq 180 °$ gradually varies, undergoes bifurcation behavior when K_b is fixed. Consequently, a set of simulations is conducted to demonstrate the changes in the $ξ$ versus $δ$ bifurcation diagrams as the joint angle $δ$ and K_b gradually varies at a consistent input speed of 62 rad/s.

Figure 6.

A series of hinge bifurcation diagrams when: (a) $δ$ = 20°, (b) $δ$ = 90°, (c) $δ$ = 160°, and (d) $δ$ = 180°.

Figure 6 represents a series of graphs showing the bifurcation behavior of the arm joint system at different angles of misalignment. Bifurcation diagrams are used to visualize changes in the dynamics of a system when a parameter, such as a misaligned angle, varies. From Figure 6, it is evident that various bifurcation types occur before the onset of chaos. A pitchfork bifurcation is observed in Figure 6(a) $δ$ = 20°, the unstable state follows the lower branch until reaching $ϕ_{1}$ = 3 rad, leading to an abrupt transition to the alternative steady state and maintaining a steady state at $ξ$ = 1 rad. There is an approximate significant region of bistability for 0 < $ϕ_{1}$ < 3 approximately in Figure 6(a), where the symmetric behavior shifts into an asymmetric oscillatory state. Subsequently, in Figure 6(b) at $ϕ_{1}$ = 5.2 rad, within the range of 20° < $δ$ < 90°, a period-doubling transition occurs, causing the period-2 orbit to bifurcate into a period-4 orbit. This period-doubling cascade continues, resulting in a broad chaotic window. This window then abruptly bifurcates into a period-2 orbit in Figure 6(c) at $ϕ_{1}$ = 6.5 rad, when 90° < $δ$ < 160°, initiating another set of period-doubling sequences that expand the chaotic domain for ϕ₁ =7–7.8 rad. It is noteworthy that in Figure 6(d), for $ϕ_{1}$ = 8.5–9.8 rad when 160° < $δ$ < 180°, the system exhibits transient chaos characterized by a sort of symmetrical periodic nature, a characteristic feature of the hard Duffing equation, as reported by Pezeshki and Dowell.²⁷ The chaotic domain exhibits a period-doubling sequence, resulting in harmonic oscillations in which the cyclic behavior is repeated consistently, as shown in Figure 6(d).

Limitations of the perturbation model

Based on the series of bifurcation diagrams, when $δ$ > 120°, the perturbation function exhibits a diverse range of behaviors. The perturbation model’s reliability is compromised in the range of $160^{0} < δ < 180^{0}$ as the stronger hands may not endure the twisting motion and the hysteresis cycle tends to repeat. However, it remains valuable for illustrating intricate growth patterns arising from basic assumptions.

Observation 2: This section discusses the sensitivity of a system to the threshold load and its potential for chaos and bifurcation behavior based on parameter values. The dynamics of the system are highly sensitive, especially at low K_b for the DC motor and high angular deflection $δ$ , highlighting behavioral differences among hinge joint models. The perturbation function remains valid within the linear range of the arm-rod coupling, but as approaches the critical value of 120°, the system transitions from quasiperiodic to chaotic bifurcation. After the joint fluctuation, the vibration amplitude experiences a varying reduction in both chaotic and quasiperiodic motion, influenced by the added motor parameters. This observation indicates that the joint-induced disturbance contributes to some extent to partially attenuate the forced vibration of the DC motor in response to external excitation.

In the subsequent section, particular attention is directed toward the integration of the electromechanical control system with the hand-arm system under the influence of the controller parameter. Consequently, the combination of the motor, gear, and hand-arm components, illustrated in Figure 7, will generate the system’s response to the impulsive disk disturbance.

Figure 7.

(a) Schematic description of the mechanical arm system and (b) conceptual design of the mass spring damper.

Mathematical modeling of the system

System description and assumptions

The biomechanical model consists of three discrete masses, dampers, and springs (Figure 7). The model assumes that the mass of each rod is concentrated at its center of gravity and that rotational inertia is neglected. It also assumes that the upper and lower rods have significant angular misalignment due to the coupling of a single U-joint to maintain the twisting motion of the rods and preserve nonlinear coupling effects, as shown in Figure 7(a).

Following the works of Tchomeni and Alugongo,²⁸ the assumption is made that the elbow’s d.o.f joint angle $δ$ is modeled as an entire joint, represented as a universal joint. The coupler shaft is assumed to be massless and its flexural vibrations are mitigated by appropriately placed idealized torsional dampers as illustrated in Figure 7(b). Consequently, this study focuses exclusively on the analysis of torsional vibrations within the system. The d.o.f.s of the system are three angular deflections of the geometric center of the arm ( $φ_{1}$ , $φ_{2}$ , and $φ_{3}$ ), two torsional deflections of the rigid bodies ( $ϕ_{1}$ and $ϕ_{2}$ ), and one lateral deflection X. To calculate the total energy of the system, first, the center of gravity of the arm rods and hinges, with the top and bottom rods representing arms 1, 2, and 3, respectively, in the Cartesian plan is determined as:

\begin{matrix} {\vec{OG}}_{1} = {\begin{matrix} x_{G 1} (t) = X (t) - a_{1} \sin φ_{1} (t) \\ {y_{G}}_{1} (t) = - a_{1} \cos φ_{1} (t) \end{matrix} \end{matrix}

(11)

\begin{matrix} {\vec{OG}}_{2} = {\begin{matrix} x_{G 2} (t) = X (t) - L_{1} \sin φ_{1} (t) - a_{2} \sin (φ_{1} (t) - φ_{2} (t)) \\ {y_{G}}_{2} (t) = - L_{1} \cos φ_{1} (t) - a_{2} \cos (φ_{1} (t) - φ_{2} (t)) \end{matrix} \end{matrix}

(12)

The center of gravity of the mass of link-1 to 3 can be expressed as:

\begin{matrix} {\vec{OG}}_{h g_{1}} = {\begin{matrix} x_{Gh g_{1}} (t) = X (t) - L_{1} \sin φ_{1} (t) \\ {y_{G}}_{h g_{1}} (t) = - L_{1} \cos φ_{1} (t) \end{matrix} \end{matrix}

(13)

\begin{matrix} {\vec{OG}}_{h g_{2}} = {\begin{matrix} x_{Gh g_{2}} (t) = X (t) - L_{2} \cos φ_{2} (t) + b \cos φ_{3} (t) \\ {y_{G}}_{h g_{2}} (t) = b \sin φ_{3} (t) \end{matrix} \end{matrix}

(14)

The linear velocity of each component can be determined by taking derivatives of equations (11)–(14) over time. However, deriving the governing equation for a coupled double-arm system is challenging due to intricate derivative calculations. The Lagrangian energy approach, which analyzes the kinetic, potential, and Rayleigh dissipation energies, offers a more efficient method to formulate the system equations.

Energy principle expression

The sum of the kinetic energies of translation and rotation of the geared motor is given by the following.

T_{MG} = \frac{1}{2} (M_{m} + \frac{J_{m} N_{g}^{2}}{r_{mg}^{2}}) \overset{\cdot}{X} (t)^{2}

(15)

The combined rotational inertia of the gears at the input link is grouped and represented by J_m, and r_mg is the load on the output link simulated through viscous torsional resistance. For real-time implementation, it is assumed that the length L₂ of the output shaft is slightly longer than that of the input shaft L₁, both being treated as rigid shafts. The total kinetic energy due to motion is the sum of the translational and rotational energies of each component, and it is expressed as:

\begin{matrix} T_{Total} = \frac{1}{2} (J_{b 1} + J_{D} + M_{1} a_{1}^{2} + M_{2} a_{2}^{2} + M_{hg} L_{1}^{2} + M_{2} L_{1}^{2} + 2 M_{2} L_{1} a_{2} \cos (φ_{2} (t))) {\overset{\cdot}{φ}}_{1} (t)^{2} \\ + \frac{1}{2} (M_{1} + M_{2} + M_{hg} + M_{hg 2} + M_{m} + \frac{J_{m} N_{g}^{2}}{r_{mg}^{2}}) \overset{\cdot}{X} (t)^{2} + \frac{1}{2} (J_{b 1} + J_{b 2} Γ {(t)}^{2}) {\overset{\cdot}{ϕ}}_{1} (t)^{2} \\ + \frac{1}{2} (J_{b 2} + J_{D} + M_{2} a_{2}^{2} + M_{hg 2} L_{2}^{2} \sin (φ_{2} (t))) {\overset{\cdot}{φ}}_{2} (t)^{2} + \frac{M_{hg 2} b^{2}}{2} {\overset{\cdot}{φ}}_{3} (t)^{2} \\ - M_{2} a_{2} (a_{2} + L_{1} \cos (φ_{2} (t))) {\overset{\cdot}{φ}}_{2} (t) {\overset{\cdot}{φ}}_{1} (t) + (M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) + M_{hg 2} L_{2} \sin (φ_{2} (t))) {\overset{\cdot}{φ}}_{2} (t) \overset{\cdot}{X} (t) \end{matrix}

\begin{matrix} - [M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) + (M_{2} L_{1} + M_{hg} L_{1} + M_{1} a_{1}) \cos φ_{1} (t)] {\overset{\cdot}{φ}}_{1} (t) \overset{\cdot}{X} (t) \\ - M_{hg 2} b \sin (φ_{3} (t)) {\overset{\cdot}{φ}}_{3} (t) \overset{\cdot}{X} (t) - M_{hg 2} L_{2} b \sin (φ_{3} (t)) \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{3} (t) {\overset{\cdot}{φ}}_{2} (t) \end{matrix}

(16)

where J_b₁, J_b₂, and J_D are, respectively, the moment of inertia of the first and second arm rods and the moment of inertia of the concentrated mass at their pivot point O.

J_{b 1} = \frac{1}{12} M_{1} {(L_{1} + a_{1})}^{2}, J_{b 2} = \frac{1}{12} M_{2} (L_{1} + L_{2})^{2} + M_{2} (L_{1} + a_{2})^{2}, J_{D} = \frac{1}{2} M_{D} (L_{1} + L_{2} + b)^{2}

(17)

The hinge between shaft 1 and shaft 2 is modeled as a twisting motion, where $k_{t 1} = k_{t 2}$ and $k_{τ 1} = k_{τ 2}$ represents the torsional spring constants of the shaft between the components. The potential energies of the rod arms and flexible hinges associated with the torsional spring stiffness and the total potential energy of the system are given below:

\begin{matrix} V_{T} (X, φ_{1}, φ_{2}, ϕ_{1}, φ_{3}) = \frac{1}{2} k_{t 2} {(φ_{2} (t) - φ_{1} (t))}^{2} + \frac{1}{2} k_{τ 1} {(ε (t))}^{2} + \frac{1}{2} k_{τ 2} {(φ_{3} (t) - φ_{2} (t))}^{2} \\ + \frac{1}{2} k_{t 1} X (t)^{2} + M_{2} g a_{2} \cos (φ_{2} (t) - φ_{1} (t)) + (M_{1} g a_{1} + M_{hg} g L_{1} + M_{2} g L_{1}) \cos φ_{1} (t) \end{matrix}

(18)

The corresponding Rayleigh dissipation function used for this system is given as:

R_{T} = \frac{1}{2} c_{t 1} \overset{\cdot}{X} (t)^{2} + \frac{1}{2} c_{τ 1} {[(1 - Γ (t)) {\overset{\cdot}{ϕ}}_{1} (t)]}^{2} + \frac{1}{2} c_{t 2} {({\overset{\cdot}{φ}}_{2} (t) - {\overset{\cdot}{φ}}_{1} (t))}^{2} + \frac{1}{2} c_{τ 2} {({\overset{\cdot}{φ}}_{3} (t) - {\overset{\cdot}{φ}}_{2} (t))}^{2}

(19)

Thus, the torsional damping coefficients of the hinge of the bonding coupler are $c_{t 1} = c_{t 2}$ and $c_{τ 1} = c_{τ 2}$ . Taking the derivatives of equation (16) through equation (19) with respect to the generalized coordinate X and rearranging the movement of the arms follows a nonlinear behavior modeled as follows:

\begin{matrix} (M_{1} + M_{2} + M_{hg} + M_{hg 2} + M_{m} + \frac{J_{m} N_{g}^{2}}{r_{mg}^{2}}) \overset{\cdot\cdot}{X} (t) \\ - ((M_{1} a_{1} + M_{2} L_{1} + M_{hg} L_{1}) \cos φ_{1} (t) + M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t))) {\overset{\cdot\cdot}{φ}}_{1} (t) \\ + (M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) + M_{hg 2} L_{2} \sin (φ_{2} (t))) {\overset{\cdot\cdot}{φ}}_{2} (t) \\ - M_{hg 2} b \sin φ_{3} (t) {\overset{\cdot\cdot}{φ}}_{3} (t) + c_{t 1} \overset{\cdot}{X} + k_{t 1} X = V_{cX} (t) \end{matrix}

(20)

Similarly, the derivatives of equations (16)–(19) with respect to $φ_{1}$ , $ϕ$ , $φ_{2}$ , $φ_{3}$ respectively yields:

\begin{matrix} - ((M_{2} L_{1} + M_{hg} L_{1} + M_{1} a_{1}) \cos φ_{1} (t) + M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) \overset{\cdot\cdot}{X} (t) \\ + (J_{D} + J_{b 1} + M_{1} a_{1}^{2} + M_{2} a_{2}^{2} + M_{2} L_{1}^{2} + M_{hg} L_{1}^{2} + 2 L_{1} M_{2} a_{2} \cos (φ_{2} (t))) {\overset{\cdot\cdot}{φ}}_{1} (t) \\ - M_{2} a_{2} (a_{2} + L_{1} \cos (φ_{2} (t))) {\overset{\cdot\cdot}{φ}}_{2} (t) + c_{t 2} ({\overset{\cdot}{φ}}_{1} (t) - {\overset{\cdot}{φ}}_{2} (t)) + k_{t 2} (φ_{1} (t) - φ_{2} (t)) = V_{c φ_{1}} (t) \end{matrix}

(21)

(J_{b 1} + J_{b 2} Γ_{1} {(t)}^{2}) {\overset{\cdot\cdot}{ϕ}}_{1} (t) + \frac{c_{τ 1} Γ_{3} (t)}{Γ_{2} {(t)}^{2}} {\overset{\cdot}{ϕ}}_{1} (t) - \frac{k_{τ 1} Γ_{3} (t) ε (t)}{Γ_{2} (t)} ϕ_{1} (t) = V_{c ϕ_{1}} (t)

(22)

\begin{matrix} (M_{hg 2} L_{2} \sin (φ_{2} (t)) + M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t))) \overset{\cdot\cdot}{X} (t) \\ - (M_{2} a_{2} (a_{2} + L_{1} \cos (φ_{2} (t))) {\overset{\cdot\cdot}{φ}}_{1} (t) + (J_{D} + J_{b 2} + M_{2} a_{2}^{2} + M_{hg 2} L_{2}^{2} \sin {(φ_{2} (t))}^{2}) {\overset{\cdot\cdot}{φ}}_{2} (t) \\ - M_{hg 2} L_{2} b \sin φ_{3} (t) \sin (φ_{2} (t)) {\overset{\cdot\cdot}{φ}}_{3} (t) + c_{t 2} ({\overset{\cdot}{φ}}_{2} (t) - {\overset{\cdot}{φ}}_{1} (t)) + c_{τ 2} ({\overset{\cdot}{φ}}_{2} (t) - {\overset{\cdot}{φ}}_{3} (t)) \\ + k_{t 2} (φ_{2} (t) - φ_{1} (t)) + k_{τ 2} (φ_{2} (t) - φ_{3} (t)) = V_{c φ_{2}} (t) \end{matrix}

(23)

\begin{matrix} - M_{hg 2} b \sin (φ_{3} (t)) \overset{\cdot\cdot}{X} (t) - M_{hg 2} L_{2} b \sin (φ_{3} (t)) \sin (φ_{2} (t)) {\overset{\cdot\cdot}{φ}}_{2} (t) + M_{hg 2} b^{2} {\overset{\cdot\cdot}{φ}}_{3} (t) \\ + c_{τ 2} ({\overset{\cdot}{φ}}_{3} (t) - {\overset{\cdot}{φ}}_{2} (t)) + k_{τ 2} (φ_{3} (t) - φ_{2} (t)) = V_{c φ_{3}} (t) \end{matrix}

(24)

Thus, the combined equations (20)–(24) are the equations of motion of a coupled arm system, established using the principles of engineering dynamics, relevant and suitable for the analysis of the behavior and performance of the system in practical engineering applications. However, it can be noted that equation (24) demonstrates that vibration in the $ϕ$ direction is dissociated from vibration in other coordinates since their equations do not share common variables for absorption. The governing equations of motion, which describe the behavior of the entire coupled arm model, can be expressed in matrix form for enhanced clarity and ease of manipulation as:

\begin{matrix} [\begin{matrix} m_{XX} & m_{X φ_{1}} & 0 & m_{X φ_{2}} & m_{X φ_{3}} \\ m_{φ_{1} X} & m_{φ_{1} φ_{1}} & 0 & m_{φ_{1} φ_{2}} & 0 \\ 0 & 0 & J_{b 1} + J_{b 2} Γ^{2} & 0 & 0 \\ m_{φ_{2} X} & m_{φ_{2} φ_{1}} & 0 & m_{φ_{2} φ_{2}} & m_{φ_{2} φ_{3}} \\ m_{φ_{3} X} & 0 & 0 & m_{φ_{3} φ_{2}} & m_{φ_{3} φ_{3}} \end{matrix}] {\begin{matrix} \overset{\cdot\cdot}{X} \\ {\overset{\cdot\cdot}{φ}}_{1} \\ \overset{\cdot\cdot}{ϕ} \\ {\overset{\cdot\cdot}{φ}}_{2} \\ {\overset{\cdot\cdot}{φ}}_{3} \end{matrix}} + [\begin{matrix} c_{t 1} + Δ_{mg} & 0 & 0 & 0 & 0 \\ 0 & c_{t 2} & 0 & - c_{t 2} & 0 \\ 0 & 0 & c_{τ 1} & 0 & 0 \\ 0 & - c_{t 2} & 0 & c_{t 2} + c_{τ 2} & - c_{τ 2} \\ 0 & 0 & 0 & - c_{τ 2} & c_{τ 2} \end{matrix}] {\begin{matrix} \overset{\cdot}{X} \\ {\overset{\cdot}{φ}}_{1} \\ \overset{\cdot}{ϕ} \\ {\overset{\cdot}{φ}}_{2} \\ {\overset{\cdot}{φ}}_{3} \end{matrix}} \\ + [\begin{matrix} k_{t 1} & 0 & 0 & 0 & 0 \\ 0 & k_{t 2} & 0 & - k_{t 2} & 0 \\ 0 & 0 & k_{τ 1} ζ (t) & 0 & 0 \\ 0 & - k_{t 2} & 0 & k_{t 2} + k_{τ 2} & - k_{τ 2} \\ 0 & 0 & 0 & - k_{τ 2} & k_{τ 2} \end{matrix}] {\begin{matrix} X \\ φ_{1} \\ ϕ \\ φ_{2} \\ φ_{3} \end{matrix}} + {\begin{matrix} V_{cX} \\ V_{c φ_{1}} \\ V_{c ϕ} \\ V_{c φ_{2}} \\ V_{c φ_{3}} \end{matrix}} = {\begin{matrix} R_{X} \\ 0 \\ 0 \\ 0 \\ 0 \end{matrix}} \end{matrix}

(25)

\begin{matrix} m_{XX} & = M_{1} + M_{2} + M_{hg} + M_{hg 2} + M_{m} + J_{m} N_{g}^{2} / r_{mg}^{2}; \\ m_{X φ_{1}} & = m_{φ_{1} X} = - (M_{1} a_{1} + M_{2} L_{1} + M_{hg} L_{1}) \cos φ_{1} (t) - M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) \\ m_{X φ_{2}} & = m_{φ_{2} X} = M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) + M_{hg 2} L_{2} \sin (φ_{2} (t)); m_{X φ_{3}} = m_{φ_{3} X} = - M_{hg 2} b \sin (φ_{3} (t)) \\ m_{φ_{1} X} & = - (M_{2} L_{1} + M_{hg} L_{1} + M_{1} a_{1}) \cos (φ_{1} (t)) - M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)) \end{matrix}

(26)

\begin{matrix} m_{φ_{1} φ_{1}} & = J_{D} + J_{b 1} + M_{1} a_{1}^{2} + M_{2} a_{2}^{2} + M_{2} L_{1}^{2} + M_{hg} L_{1}^{2} + 2 L_{1} M_{2} a_{2} \cos (φ_{2} (t)); \\ m_{φ_{2} X} & = M_{hg 2} L_{2} \sin (φ_{2} (t)) + M_{2} a_{2} \cos (φ_{2} (t) - φ_{1} (t)); m_{φ_{1} φ_{2}} = m_{φ_{2} φ_{1}} = - M_{2} a_{2} (a_{2} + L_{1} \cos (φ_{2} (t))); \\ m_{φ_{2} φ_{2}} & = J_{D} + J_{b 2} + M_{2} a_{2}^{2} + M_{hg 2} L_{2}^{2} \sin (φ_{2} (t))^{2}; m_{φ_{2} φ_{3}} = m_{φ_{3} φ_{2}} = - M_{hg 2} L_{2} b \sin (φ_{3} (t)) \sin (φ_{2} (t)) \\ m_{φ_{3} φ_{3}} & = M_{hg 2} b^{2}; and Δ_{mg} = K_{t} N_{g} K_{b} / R_{a} r_{mg}^{2}; Γ_{d} = Γ_{3} / {Γ_{2}}^{2}; Γ_{d 1} = - Γ_{3} / Γ_{2} \end{matrix}

(27)

Γ_{2} = N^{2} (γ^{2} N^{2} + γ^{2} + 4 γ + 1) + 1; Γ_{3} = (γ + 1) (γ N^{2} - 1) (N^{2} + 1), N = \tan ϕ_{1} (t), γ = \cos δ

(28)

R _X and Vc_n represent the generalized vector forces and some of the nonlinear terms derived from the coupling generalized coordinates. The specific expressions for each of these generalized vector terms are provided in Appendix A.

Numerical results of dynamical non-linear systems

As an initial objective of the overall analysis of the device, it is necessary to solve the equation that governs the oscillation of the system. Equation (25) is used to study system stability control and harmonic solutions using fast Fourier transform, orbit, and arm displacement, thus providing real-time feedback on weight variation. The study involves numerical solutions for a nonlinear model representing a coupled arm system, uncovering chaotic fluctuations under specific conditions using Poincaré maps and bifurcation diagrams. A representative set of parameter values described in Table 2 is chosen to perform simulations through the time domain expressions of equation (25).

Table 2.

Coupled arms generating system parameters.

Symbol	Description	Value and units
c _τ1	Torsional viscous damping coefficient of the hinge 1	24 N.m.s/rad
c _τ2	Torsional viscous damping coefficient of the Hinge 2	5.5 N.m.s/rad
c _t1	Viscous damping coefficient of the upper arm 1	2.0 × 10³ ns/m
c _t2	Viscous damping coefficient of the lower arm 2	138.42 ns/m
g	Gravitational constant	9.81 m/s²
J _b1	Moment of inertia of the upper arm rod at its center of gravity	26.347 kg/m²
J _b2	Moment of inertia of the lower arm rod at its center of gravity	24.05 kg/m²
J _m	Rotational moment of inertia of the DC Motor’s input shaft	3.9 × 10⁻³ kg/m
J _D	Moment of inertia of the external disk weight mass	24.05 kg/m²
k _τ1	Upper arm 1 torsional stiffness	2.5 × 10⁵ N/m.rad
k _τ2	Lower arm 2 torsional stiffness	2.5 × 10⁵ N/m.rad
K _t1	Upper arm 1 stiffness	6 × 10⁷ N/m
K _t2	Lower arm 2 stiffness	6 × 10⁷ N/m
a ₁	Length of upper arm 1 from pivot 1 to center of gravity	0.1143 m
a ₂	Length of lower arm 2 from pivot 2 to center of gravity	0.1778 m
L ₁	Length of upper arm 1	0.209 m
L ₂	Length of lower arm 2	0.3365 m
b	Length of hand	0.2165 m
r _mg	Motor pinion radius	0.063 m
M _m	Motor mass	10.57 kg
M _hg1	Hinge mass L₁ to L₂	1.178 kg
M _hg2	Hinge mass L₂ to b	0.178 kg
M _D	External disk weight mass	Adjustable kg
M ₁	Upper arm 1 mass	0.072 kg
M ₂	Lower arm 2 mass	0.127 kg
Ra	Motor armature resistance	2.6 Ω
Ng	Planetary gearbox gear ratio	3.71

Numerical results and discussion

Due to the strongly nonlinear characteristics of the joint impact arm system, extracting precise vibration responses directly poses a challenge. Obtaining exact solutions proves to be difficult, leading to the transformation of the derived equation (25) of motion into a set of first-order differential equations. Subsequently, the fourth-order Runge-Kutta method is applied to solve these equations and acquire data. In the course of the calculation, a smaller integration step is selected to guarantee the stability of the solution and prevent numerical divergence, particularly at points where the derivatives of the impact forces of the load are present.

(a) First case study: Study of motion for M_D = 0, K_b = 1.982, $δ$ = 60°, and to accommodate movement or dampen vibration, especially under load conditions, (kτ₁ >> K_t1 and kτ₂ >> K_t2). The dynamic response of the electromechanical system is obtained by solving equation (25) using Matlab. Figure 8 illustrates the frequency evolution of $φ_{1}$ , $φ_{2}$ , and $φ_{3}$ when no external weight is applied to the coupled link system. The vibration signals in the frequency domain are acquired by applying the fast Fourier transform (FFT) on the signals in the time domain. Analysis of Figures 8 –11 indicates that each mechanical arm frequency moves with a different amplitude.

Figure 8.

Amplitude frequency response for each link arm at a critical point (M_D = 0, K_b = 1.982, $δ$ = 60°): (a) frequency spectrum, (b) 3D waterfall transform, (c) Poincaré maps for system, and (d) projection of the Poincaré maps expansion.

Figure 9.

Amplitude frequency response for each link arm at a critical point (M_D = 15, K_b = 1.982, $δ$ = 60°): (a) frequency spectrum, (b) 3D waterfall transform, (c) Poincaré maps for system, and (d) projection of the Poincaré maps expansion.

Figure 10.

Amplitude frequency response for each link arm at the critical threshold (M_D = 35 kg, K_b = 1.982, $δ$ = 60°): (a) frequency spectrum, (b) 3D waterfall transform, (c) Poincaré maps for system, and (d) projection of the Poincaré maps expansion.

Figure 11.

Amplitude frequency response for each arm above the critical threshold (M_D = 50 kg, K_b = 1.982, $δ$ = 60°): (a) frequency spectrum, (b) 3D waterfall transform, (c) Poincaré maps for system, and (d) projection of the Poincaré maps expansion.

Aperiodic (chaotic) behavior in a dynamical system is characterized by a broadband frequency spectrum. In addition, the frequency spectrum is densely populated, with multiple peaks observed around the fundamental frequency. The electromechanical system exhibits significant resonance peaks near the first critical speed. Observations of the frequency spectrum in Figures 8(a), 9(a), 10(a), and 11(a) reveal peaks at frequencies 42.123, 130.381, and 764.232 Hz corresponding to spindle speed rotation. Figures illustrate that arm 3 is significantly more affected by the preload action compared to input arm 1 and intermediate arm 2. Figure 8(a) shows that the third arm of the system displays super-harmonic resonances (also referred to as subcritical resonances) equivalent to the combination of the input link 1 and intermediate link 2’s first-order critical speed. Under a charging load of 50 kg, high amplitudes were observed at the same resonant frequency points. However, monitoring the evolution of the frequencies for each arm does not indicate a malfunction in the system, and there is no evidence of fault frequencies for the analyzed arms. It implies that a system with three fundamental frequencies is referred to as a third-order periodic system.

(b) Second case study: Study of motion for the control parameter M_D = 15 kg, K_b = 1.982, $δ$ = 60°.

The spectrum presented in Figure 8(a) lacks a rich variable content frequency. The observed divergence in frequency amplitudes indicates the presence of the jump phenomenon called “resonance.” To illustrate nonlinear resonance, temporal simulations were impractical due to system sensitivity. The 3D waterfall power spectrum technique is used to identify various types of chaotic behavior in systems. For nonlinear systems, a multiple of fundamental frequencies indicates superharmonics, while incommensurable ratios suggest quasiperiodic behavior. The aperiodic behavior is characterized by broadband frequency spectra and dense frequency spectra, with more peaks around the fundamental frequency. Instead, a 3D waterfall transform analysis revealed systematic variations in frequency over time. Figures 8(b), 9(b), 10(b), and 11(b) display the 3D waterfall wavelet synchrosqueezing distribution of the signals from the three arms at different speeds. In contrast to the previous FFT spectrum, it is now much easier to identify the type of chaotic motion visible in the signal. The frequency-based disturbance analysis can be distinguished among other frequencies in this case. The 3D waterfall representation allows for clear identification of time and frequency for each frequency, making it easy to interpret and confirm frequency complexity (number of harmonics and increases in amplitude at rotation frequencies), indicating system instability. When a very light mass is introduced, a relatively rapid increase in the period is observed. For M_D = 15 kg, the motion, although nonchaotic, tends to become increasingly unpredictable. This trend toward chaos is confirmed with a Poincaré section where the initial conditions are slightly different to observe the presence or absence of a major characteristic of chaos. Sensitivity to the threshold load is obtained by modifying the value of the mass to be transported by arm 3 (M_D = 35 and 50 kg).

Poincaré maps are used to simplify complex systems and provide insight into their nonlinear dynamic behavior. A closed curve on the map indicates quasiperiodic behavior, making it possible to visually distinguish nonchaotic movements. The initiation of chaos on the Poincaré map is marked by a distinct set of points. For M_D = 35 kg, a distinctive phenomenon of recurrence of nonsimple periodicity is observed. The Poincaré section for this value exhibits an ellipsoidal curve, as indicated by the first Figures 8(c), 9(c), 10(c), and 11(c). The subsequent image, depicting the Poincaré section for M_D = 35 kg, illustrates a less predictable behavior. However, the overlay of the two solutions with different threshold loads eliminates the possibility of chaos. Once again, the sensitivity of the movement is influenced by various values of the M_D parameter that can be examined. In this context, different behaviors are observed more broadly. It is evident that around M_D = 35 kg, the movement becomes more responsive to the threshold load, progressively becoming more unstable until M_D = 50 kg, eventually transitioning into chaos. The region near M_D = 50 kg signifies a shift toward chaos, as solutions with close initial conditions diverge, leading to a noticeable disparity between scattered points (Figures 8(d), 9(d), 10(d), and 11(d). It is observed that a high value of M_D significantly influences the movement, leading it toward greater chaos. This intuition is supported by a Poincaré and bifurcation section, specifically for M_D = 50 kg while maintaining the angle and K_b constant. The observation here indicates that, for a substantial M_D, in contrast to the previous scenario, the movement becomes highly chaotic, with the chaos manifesting definitively and prominently.

(d) Fourth case study: Study of motion for the control parameter M_D = 50 kg, K_b = 1.982, $δ$ = 60°.

Observation 3: Quasi-periodic motion exhibits high peaks at the natural frequency and harmonics, while chaotic motion displays a wide frequency range and disturbances. Figures 8 –11 reveal a diverse range of behaviors in the Poincaré section as the motion begins, aiding in the clarification of complex trajectories within higher-dimensional spaces. As depicted in Figures 8(c), 9(c), 10(c), and 11(c), as well as Figures 8(d), 9(d), 10(d), and 11(d), the evolving Poincaré section illustrates the transformation as the oscillation level increases. There is a clear distinction between periodic behavior and quasiperiodic behavior, characterized by a very long period cycle that never closes upon itself. In Figures 8(c), 9(c), 10(c), and 11(c), the quasiclosed orbits, indicative of quasiperiodic behavior, gradually evolve into irregular patterns as the number of iterations systematically tends toward infinity. The distinction between chaotic and other movements is visually perceptible on the Poincaré map, quasi-periodic behavior appears as closed curves or points, and the presence of a distinct set of points indicates chaos.

Bifurcations of nonlinear arm systems for different mass values

To conclude, four bifurcation diagrams will be generated, each corresponding to different values of the external weight load, M_D. According to the specific magnitudes of the M_D external weight load, the bifurcation analysis of the nonlinear system is performed and presented in Figure 12.

Figure 12.

Bifurcation diagram of equation (25) when K_b = 1.982, $δ$ = 60° and changes in the adjustable M_D: (a) parameter K_b = 1.982 and M_b = 0 kg, (b) parameter K_b = 1.982 and M_b = 15 kg, (c) parameter K_b = 1.982 and M_b = 0 kg, (c) parameter K_b = 1.982 and M_b = 35 kg, and (d) parameter K_b = 1.982 and M_b = 50 kg.

During this examination, the parameter K_b is kept constant at a value of K_b = 1.982, while and = 60°. The values chosen for M_D in this section are M_D = 0, 15, 35, and 50 kg, allowing a thorough analysis of the behavior of the system in a range of external weight loads. The chaotic behavior of the system is also validated by the bifurcation diagram, shown in Figure 12, where the chaotic behavior is based on multiple irregularities, seemingly random frequencies, and a lack of clear patterns that can correspond to different modes or components of the coupled arms. In system (25), the amplitude of the control input, M_D, is regarded as a bifurcation parameter, while the others are kept fixed, and it is calculated to illustrate the various dynamical behaviors of equation (25). The bifurcation diagram for parameter M_D reveals both chaotic and periodic behavior of the system across a broad range of M_D values, exhibiting a stable limit cycle followed by transient chaos. Using these four bifurcation diagrams, it is observed that the impact of forced motion on the chaotic nature of movement depends on its relative intensity with respect to M_D and the inherent stability of motion without external influence. Specifically, it is observed that for 0 < M_D < 50 kg, the system exhibits a highly periodic motion in the absence of forced movement, and chaos requires M_D to exceed 35 kg. This is comprehensible, given our earlier conclusion that, for forced movement, M_D > 35 kg induces system unbalance, establishing 35 kg as the threshold value for chaos.

Observation 4: Various tools have been employed to analyze the complex dynamic behavior of the proposed coupled arms manipulator, such as the Poincaré map, the bifurcation diagram, and the frequency spectrum. However, bifurcation and Poincaré analyses offer a more comprehensive understanding of system fluctuations, especially regarding stability for specific joint angles. Unlike classical techniques, which may not capture the nuanced behavior of the system, these analyses demonstrate pronounced link stability until a certain mass-amplitude threshold is exceeded. The system exhibits chaotic behavior with a stable limit cycle in its steady state, which is influenced by changes in the amplitude of the control input. The bifurcation analysis of the system is characterized by the existence of a folding, which means a single equilibrium zone, while elsewhere there is no equilibrium point. In Figure 12, the black points correspond to stable zones, and the red and green points to unstable zones. Taking into consideration the value of the external load weight, the stable zone is considerably reduced and eventually widens so that there is a large zone of divergent points. This oscillation zone is unstable: The loaded arm system evolves into a chaotic zone. The chaotic nature of the proposed system is further validated in the figure below through the examination of time displacement and orbit patterns for each arm, revealing both single- and multi-circle chaotic orbits.

Dynamic analysis based on arms response

Considering the topic of this article, one of the objectives is to provide a qualitative interpretation of the results and an analysis of the behavior of interconnected arms from the perspective of bifurcation theory. In dedicated terminology, the problem has external load M_D as a control parameter and gain K_b = 1.982 as a fixed parameter corresponding to a bifurcation of nodes of periodic orbits. It remains to confirm this assertion through two-time simulations initialized with values much larger than the amplitude of the stable periodic orbit so that they converge toward it if it exists. In this section, the nonlinear dynamics of arms connected via a flexible joint are analyzed, revealing complex motion patterns, oscillation behaviors, and stability characteristics under various loads and system configurations. A control gain K_b = 1.982 is selected and the simulation is initialized with a minimal lifting mass M_D = 0 kg and a maximum mass. M_D < 50 kg.

Observation 5: Parametric studies were carried out to study the evolution of the bifurcation points and the influence of different parameters on the instability zone. Orbit observations have contributed to the formulation of simulation strategies that favor the observation of limit cycles and their characteristics. However, at high loads, it is impossible to capture a periodic orbit because the time simulations never actually pass through the same point. It is therefore difficult not only to draw conclusions about the real existence of periodic orbits but also to have a framework for studying their characteristics. The observation of the diagrams presented in Figure 13(a) and (b) allows for the estimation of a stable limit cycle, regardless of the initial disturbance. The influence of the variable M_D is shown in Figure 13(c) and (d) is that the amplitude of oscillations and disturbance increases as the external load increases. However, it is difficult to provide a precise estimate of their movement amplitude in order to characterize a bifurcation. For future work, the performance of experimental analyzes will provide empirical evidence to strengthen the conclusions drawn from the numerical simulations. Furthermore, it is necessary to explore and extend Lyapunov’s theorems in future contributions, especially in contexts where control dynamics will be incorporated.

Figure 13.

Time simulation of connected arms lifting a specific load (M_D = 0 and 50 kg and K_b = 1.982).

Conclusions

This article explores the steady-state operation of a DC motor connected to a flexible joint manipulator with a double arm, focusing on chaotic behaviors. The dynamic equations are formulated using the energy principle, and the model incorporates non-linearities during arm flexion and extension for enhanced accuracy. The motor introduces a d.o.f to simulate the clavicle’s role in arm movement, fostering interaction between two manipulators. The study demonstrates that due to its highly non-linear dynamic characteristics, a two-arm system designed for weightlifting exhibits complex phenomena, including higher subharmonic oscillations, bifurcations, and chaotic movements. The arms exhibit structural and kinematic complexities, such as offset attachments and variable link cross-sections. The system shows increased instability when external disturbances occur due to weight lifting. The chaotic nature of the system is confirmed using the Poincaré map and bifurcation diagram, highlighting the coexistence of chaotic orbits in all instances and providing a quantitative characterization of the system responses. The Poincaré section plays a crucial role in observing bifurcations and chaos of nonlinear oscillators. In the case of a heavy system, the response includes multiple high-order harmonic motions and chaos with significant contamination, which poses problems for easy distinction in both orbit plots and frequency spectra, and strengthens the conclusions drawn from the numerical simulations.

In future research, the authors aim to expand upon current research in the following areas:

Conduct experimental studies to validate the proposed biomechanical arm model under real-world conditions.

Employ advanced analytical techniques such as Lyapunov algorithms to enhance stability analysis and refine adaptive control methods, allowing the biomechanical arm system to adapt to varying loads more effectively.

Create a more faithful mechanical imitation of human musculoskeletal structures, advanced materials, and actuators with properties closer to those of biological tissues.

Footnotes

Appendix

Finally, the components of the partial derivatives of equation (25) are, given below:

(A-1)

R_{X} (t) = K_{t} K_{g} \frac{(V_{b} (t) r_{mg} - K_{b} \overset{\cdot}{X} (t))}{R_{a} r_{mg}^{2}}

(A-2)

\begin{array}{l} V_{c X} = ((M_{2} L_{1} + M_{2} a_{2} + L_{1} M_{h g}) \sin (φ_{1} (t)) - M_{2} a_{2} \sin (φ_{2} (t) - φ_{1} (t))) {\dot{φ}}_{1} {(t)}^{2} \\ - M_{h g 2} b \cos (φ_{3} (t)) {\dot{φ}}_{3} {(t)}^{2} + (M_{h g 2} L_{2} \cos (φ_{2} (t)) - M_{2} a_{2} \sin (φ_{2} (t) - φ_{1} (t))) {\dot{φ}}_{2} {(t)}^{2} \\ + 2 M_{2} a_{2} \sin (φ_{2} (t) - φ_{1} (t)) {\dot{φ}}_{2} (t) {\dot{φ}}_{1} (t) \end{array}

(A-3)

\begin{matrix} V_{c φ_{1}} = - 2 M_{2} L_{1} a_{2} \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{2} (t) {\overset{\cdot}{φ}}_{1} (t) + M_{2} L_{1} a_{2} \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{2} (t)^{2} \\ + M_{2} a_{2} \sin (φ_{2} (t) - φ_{1} (t)) {\overset{\cdot}{φ}}_{2} (t) \overset{\cdot}{X} (t) \\ + ((M_{2} L_{1} + M_{hg} L_{1} + M_{1} a_{1}) \sin (φ_{1} (t)) - M_{2} a_{2} \sin (φ_{2} (t) - φ_{1} (t))) {\overset{\cdot}{φ}}_{1} (t) \overset{\cdot}{X} (t) \\ - (M_{1} L_{1} + M_{hg} L_{1} + M_{1} a_{1}) g \sin (φ_{1} (t)) + M_{2} a_{2} g \sin (φ_{2} (t) - φ_{1} (t)) \end{matrix}

(A-4)

V_{c ϕ_{1}} = \frac{4 J_{b 2} Γ^{2} (γ + 1) (2 γ N^{2} + γ - 1)}{Γ_{2}}

(A-5)

\begin{matrix} V_{c φ_{2}} = M_{2} a_{2} \sin (φ_{2} (t) - φ_{1} (t)) {\overset{\cdot}{φ}}_{1} (t) \overset{\cdot}{X} (t) + M_{2} a_{2} L_{1} \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{2} (t) {\overset{\cdot}{φ}}_{1} (t) \\ + 2 M_{hg 2} L_{2}^{2} \cos (φ_{2} (t)) \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{2} (t)^{2} - M_{hg 2} L_{2} b \cos (φ_{3} (t)) \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{3} (t)^{2} \\ + M_{hg 2} L_{2} b \cos (φ_{2} (t)) \sin (φ_{3} (t)) {\overset{\cdot}{φ}}_{3} (t) {\overset{\cdot}{φ}}_{2} (t) - M_{2} a_{2} g \sin (φ_{2} (t) - φ_{1} (t)) \end{matrix}

(A-6)

\begin{matrix} V_{c φ_{3}} = - M_{hg 2} L_{2} b \cos (φ_{2} (t)) \sin (φ_{3} (t)) {\overset{\cdot}{φ}}_{2} (t)^{2} - M_{hg 2} b \cos (φ_{3} (t)) {\overset{\cdot}{φ}}_{3} (t) \overset{\cdot}{X} (t) \\ - M_{hg 2} L_{2} b \cos (φ_{2} (t)) \sin (φ_{2} (t)) {\overset{\cdot}{φ}}_{3} (t) {\overset{\cdot}{φ}}_{2} (t) \end{matrix}

and

(A-7)

Γ (t) = \frac{γ - γ N^{2} (N^{2} - 4) + 2 N^{2} + 2}{N^{2} (γ^{2} N^{2} + γ^{2} + 4 γ + 1) + 1}

Credit authorship contribution statement

B.X. Tchomeni: Conceptualization, mal analysis. A.A. Alugongo: Validation, investigation, resources, supervision. Writing review and editing, Validation.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is based on research supported in part by the Vaal University of Technology (VUT), South Africa, through the Department of Higher Education and Training University Capacity Development Grant.

ORCID iD

Bernard Xavier Tchomeni Kouejou

Data availability statement

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study

References

Spong

Hutchinson

Vidyasagar

Robot modeling and control. John Wiley & Sons, 2020.

Al-Qahtani

Mohammed

Sunar

Dynamics and control of a robotic arm having four links. Arab J Sci Eng 2017; 42(5): 1841–1852.

Yuguang

Fan

Dynamic modeling and adaptive fuzzy sliding mode control for multi-link underwater manipulators. Ocean Eng 2019; 187: 106202.

Duprey

Naaim

Moissenet

, et al. Kinematic models of the upper limb joints for multibody kinematics optimisation: an overview. J Biomech 2017; 62: 87–94.

Quental

Folgado

Ambrósio

, et al. A new shoulder model with a biologically inspired glenohumeral joint. Med Eng Phys 2016; 38(9): 969–977.

Bolsterlee

Veeger

HEJ

van der Helm

FCT

. Modelling clavicular and scapular kinematics: from measurement to simulation. Med Biol Eng Comput 2014; 52(3): 283–291.

Kecskeméthy

Weinberg

An improved elasto-kinematic model of the human forearm for biofidelic medical diagnosis. Multibody Syst Dyn 2005; 14(1): 1–21.

Lovasz

Mărgineanu

Ciupe

, et al. Design and control solutions for haptic elbow exoskeleton module used in space telerobotics. Mech Mach Theory 2017; 107: 384–398.

Aljinovic

Crnjac

Nikola

, et al. Integration of the human-robot system in the learning factory assembly process. Procedia Manuf 2020; 45: 158–163.

10.

Liao

Wang

Trust-based variable impedance control of human–robot cooperative manipulation. Robot Comput Integr Manuf 2024; 88: 102730.

11.

Buerkle

Eaton

Lohse

, et al. EEG based arm movement intention recognition towards enhanced safety in symbiotic human-robot collaboration. Robot Comput Integr Manuf 2021; 70: 102137.

12.

Bascetta

Ferretti

Ensuring safety in hands-on control through stability analysis of the human-robot interaction. Robot Comput Integr Manuf 2019; 57: 197–212.

13.

Faramin

Ataei

Chaotic attitude analysis of a satellite via Lyapunov exponents and its robust nonlinear control subject to disturbances and uncertainties. Nonlinear Dyn 2016; 83(1–2): 361–374.

14.

Gierlak

Warmiński

Analysis of bifurcation vibrations of an industrial robot arm system with joints compliance. Appl Sci 2023; 13(21): 11941.

15.

Madsen

Rosenlund

Brandt

, et al. Comprehensive modeling and identification of nonlinear joint dynamics for collaborative industrial robot manipulators. Control Eng Pract 2020; 101: 104462.

16.

Park

Choi

. Input-to-state stability of variable impedance control for robotic manipulator. Appl Sci 2020; 10(4): 1271.

17.

Krishan

Verma

Assessment and enhancement of hopf bifurcation stability margin in uncertain power systems. Elect Power Syst Res 2022; 206: 107783.

18.

Tartaruga

Cooper

Lowenberg

, et al. Uncertainty and sensitivity analysis of bifurcation loci characterizing nonlinear landing-gear dynamics. J Aircr 2018; 55(1): 162–172.

19.

Lee

Sifakis

Terzopoulos

. Comprehensive biomechanical modeling and simulation of the upper body. ACM T Graphic (TOG) 2009; 28(4): 1–17.

20.

Taylor

Westbrook

Head

, et al. The single-leg vertical hop provides unique asymmetry information in individuals after anterior cruciate ligament reconstruction. Clin Biomech 2020; 80: 105107.

21.

Blache

Creveaux

Dumas

, et al. Glenohumeral contact force during flat and topspin tennis forehand drives. Sports Biomech 2017; 16(1): 127–142.

22.

Safaeepour

Esteki

Ghomshe

, et al. Quantitative analysis of human ankle characteristics at different gait phases and speeds for utilizing in ankle-foot prosthetic design. Biomed Eng Online 2014; 13(1): 8.

23.

McDonald

Mulla

Keir

PJ.

Muscular and kinematic adaptations to fatiguing repetitive upper extremity work. Appl Ergon 2019; 75: 250–256.

24.

Ogata

. Modern control engineering. 5th ed. Boston: Pearson, 2010.

25.

Bernstein

. Modeling and control: applications to a double inverted pendulum and radio frequency interference. 2018.

26.

Tchomeni Kouejou

Alugongo

. Experimental analysis of a cracked Cardan shaft system under the influence of viscous hydrodynamic forces. Fluids 2023; 8(7): 211.

27.

Pezeshki

Dowell

EH.

An examination of initial condition maps for the sinusoidally excited buckled beam modeled by the Duffing’s equation. J Sound Vib 1987; 117(2): 219–232.

28.

Tchomeni