A new skeletal model for the ankle joint complex

Abstract

The talocrural and the talocalcaneal articulations collectively form the ankle joint complex of the human foot and are the focus of investigation of this work. The talocrural articulation enables plantarflexion and dorsiflexion, while the talocalcaneal articulation allows inversion and eversion of the foot. A comprehensive analysis of the literature suggests that the ankle joint complex is modeled in different manners considering approaches with varying complexity levels, which more or less accurately mimic its intrinsic anatomical features. Several studies assume that the foot articulates with the leg via the talocrural articulation only, which is modeled as a revolute joint. Other studies consider the movements allowed by both articulations and model the ankle joint complex as spherical, revolute, or classical universal joints. Most existing approaches do not consider sufficiently accurate anatomical modeling of this joint complex. Thus, this work presents a new skeletal model for the ankle joint complex of the human foot that considers the actual anatomy and movements of the talocrural and the talocalcaneal articulations. The proposed approach uses a modified universal joint, which incorporates a massless link to mimic the actual function of the talus bone. The developed formulation is compared with a model available in the literature, which uses a classical universal joint. The outcomes show that modeling the ankle joint complex as a modified universal joint allows a more realistic representation of the anatomy of the human foot. The main differences between the two joint models are observed in the mediolateral direction.

Nomenclature

Latin symbols

Symbol	Units (SI)	Description
$\mathbf{A}_{k}$	-	Rotation matrix of body k
d	m	Position vector connecting point $P_{i}$ on body i to point $P_{j}$ on body j
D	-	Jacobian matrix of the constraint equations
$e_{0}$, $e_{1}$, $e_{2}$, $e_{3}$	-	Euler parameters
g	N, N⋅m	External generalized force vector
h	-	Generic rotational joint h
i	-	Body i
I	kg⋅m2	Moment of inertia
I	-	Identity matrix
j	-	Body j
$j_{h}$	N⋅m⋅s	Damping coefficient of a generic rotational joint h
l	m	Length of a segment or link
M	kg, kg⋅m2	Mass matrix
mp,h	N⋅m	Maximum moment applied to restrict the motion of a generic rotational joint h
$\mathbf{m}_{h}^{\mathrm{d}}$	N⋅m	Joint dissipative moment vector of a generic rotational joint h
$\mathbf{m}_{h}^{\mathrm{mr}}$	N⋅m	Joint motion-restricting moment vector of a generic rotational joint h
$\mathbf{m}_{h}^{\mathrm{r}}$	N⋅m	Joint resistance moment vector of a generic rotational joint h
n	-	Unit vector normal to the surface of the circumduction cone
$O_{k}$	-	Center of mass of body k
P	-	Point representing the intersection of the axes of the classical universal joint
q	m	Vector containing the system coordinates
$\mathbf{r}_{k}$	m	Position vector of the center of mass of body k described in global coordinates
$\mathbf{r}_{k}^{P}$	m	Global position vector of point P located on body k
RoM	-	Range of motion
$\mathbf{s}_{k}^{P}$	m	Global position vector of point P located on body k with respect to local coordinates
$\mathbf{s}_{k}$	m	Vector placed on the joint axis belonging to body k
t	s	Time variable
ur,h	-	Unit vector defined using the moving body’s local reference frame
ur,h,ξ	-	Local ξ component of the unit vector ur,h
ur,h,η	-	Local η component of the unit vector ur,h
ur,h,ζ	-	Local ζ component of the unit vector ur,h
$\mathbf{u}_{\xi ,\textit{h}}$	-	Unit vector defining the ξ axis of joint h local reference frame
$\mathbf{u}_{\eta ,\textit{h}}$	-	Unit vector defining the η axis of joint h local reference frame
$\mathbf{u}_{\zeta ,\textit{h}}$	-	Unit vector defining the ζ axis of joint h local reference frame

v	m/s, rad/s	Vector containing the system velocities
$\dot{\mathbf{v}}$	m/s2, rad/s2	Vector containing the system accelerations
xyz	-	Global coordinate system
y	m, m/s, rad/s	Auxiliary vector containing the system velocities and positions
$\dot{\mathbf{y}}$	m/s, s−1,	Auxiliary vector containing the system accelerations and velocities
	m/s2, rad/s2

Greek symbols

Symbol	Units (SI)	Description
α	-	Baumgarte stabilization coefficient
β	-	Baumgarte stabilization coefficient
γ	-	Right-hand side vector of the acceleration constraint equations
ε	rad	Angle of the talocrural joint axis in the transverse plane
θ	rad	Angle of the talocalcaneal joint axis in the sagittal plane
κ	rad	Relative angular motion between the moving body and the limits of the circumduction cone
	-	Lagrange multipliers vector
	-	Right-hand side vector of the velocity constraint equations
ρ	rad	Angle of the talocalcaneal joint axis in the transverse plane
$\sigma _{h}$	rad	Latitude of a generic rotational joint h
σh,max	rad	Maximum latitude of a generic rotational joint h
$\Delta \sigma _{h}$	rad	Difference between the passive and active range of motion of a generic rotational joint h
τ	rad	Angle of the talocrural joint axis in the frontal plane
Φ	-	Position constraint equations vector
$\dot{\boldsymbol{\Phi}}$	-	Velocity constraint equations vector
$\ddot{\boldsymbol{\Phi}}$	-	Acceleration constraint equations vector
$\psi _{h}$	rad	Longitude of a generic rotational joint h
ω	rad/s	Angular velocity vector
ω̇	rad/s2	Angular acceleration vector
ξηζ	-	Body fixed coordinate system

Subscripts

Symbol	Description
0	Initial condition
h	Relative to a generic rotational joint h
i	Relative to body i
j	Relative to body j
k	Relative to body k
p	Penalty
re	Real angle
t	Relative to time instant t

Superscripts

Symbol	Description
n	Normal constraint
P	Generic point P
s	Spherical joint
u	Relative to the classical universal joint
um	Relative to the modified universal joint

Operators

Symbol	Description
( )T	Matrix or vector transpose
(′)	Components of a vector in a body-fixed coordinate system
(′′)	Components of a vector in the joint’s coordinate system
($^{\boldsymbol{\cdot}}$)	First derivative with respect to time
(⋅⋅)	Second derivative with respect to time
(∼)	Skew-symmetric matrix of a vector

1 Introduction

Over the last few years, the human movement has been a subject of extensive investigation by many authors, involving a vast interest in both clinical and sports applications [1–7]. Biomechanical models of the human body can greatly contribute to advancing knowledge in this scientific field as they provide useful information on several biomechanical parameters that are difficult or even impossible to measure in experimental settings, including joint reaction forces and moments. The validity of the predicted results strongly relies on the development of computationally efficient and anatomically accurate biomechanical models of the human system, which comprehensively model the anatomical segments, articulations, muscles, and body-environment interaction, such as the contacts that can take place between the foot and the ground surfaces [7–12]. It is important to note that the reliability of the parameters provided by biomechanical models is strongly dependent on whether these models have been validated for the intended purpose, thus preventing incorrect conclusions. Validation can be achieved through, for instance, comparing the obtained results with previously validated and published studies. Joint angles and moments-of-force, muscle activations and forces, or ground reaction forces can be utilized to validate biomechanical models [13, 14]. By providing researchers and healthcare professionals with knowledge of important biomechanical parameters, a more in-depth understanding of healthy and pathological conditions can be achieved, ultimately helping the clinical decision-making process and the development of personalized solutions specifically targeted to each patient’s needs [3, 6, 15].

Due to its intricate nature, the ankle joint complex of the human foot is quite a challenging problem to handle under the umbrella of biomechanics of motion discipline. The foot is the part of the human body that contacts the ground, and it comprises three major bone groups, namely, the tarsus, the metatarsus, and the phalanges, as can be observed in Fig. 1a [16, 17]. The tarsus group is composed of seven bones. One of them is the calcaneus, which projects posteriorly on the foot, forming the prominence of the heel. The remaining bones are the talus, the navicular, the cuboid, and the three cuneiform bones, namely medial, intermediate, and lateral cuneiforms. The metatarsus group includes five different metatarsal bones, which do not have specific names, being numbered from one to five with roman numerals, starting from the medial side of the foot. Finally, there are two phalanges in the great toe (proximal and distal), while the remaining toes comprise three phalanges each (proximal, middle, and distal) [16] (see Fig. 1a). The intricate interconnections that exist between the bones of the human foot result in the establishment of 31 articulations [17] that allow several movements during daily life activities. Amongst these 31 articulations, the talocrural and the talocalcaneal articulations, which collectively form the ankle joint complex of the human foot [18–20], are the focal points of the present research work. The talocrural articulation enables plantarflexion and dorsiflexion, and it is located between the superior surface of the talus of the foot and the distal ends of the tibia and fibula of the leg. The foot also performs inversion and eversion, which are movements provided by the talocalcaneal articulation that is located between the tarsal bones of the posterior part of the foot, more specifically between the superior surface of the calcaneus and the inferior surface of the talus [17, 21]. While the talocrural and the talocalcaneal articulations are the major contributors to plantarflexion, dorsiflexion, inversion, and eversion of the human foot, those movements result from the combination of the different degrees-of-freedom of the ankle joint complex and the foot in all cardinal planes. A schematic representation of the ankle joint complex is presented in Fig. 1b.

Fig. 1

Schematic representation of the (a) bones and (b) ankle joint complex of the right human foot (Color figure online)

The ankle joint complex of the human foot has two anatomical particularities. The first one is associated with the fact that the anatomical axes of the talocrural and the talocalcaneal articulations do not intersect at a specified point, but instead, they are separated by a certain distance. Therefore, the joint axes are non-coplanar and the distance between them can be established as the length of the talus since this bone is inserted between the two articulations (see Fig. 1b). The second particularity is that the axes of the talocrural and the talocalcaneal articulations define specific anatomical orientations between each other [22–25], as can be observed in the schematic representation of Fig. 2. Both the talocrural and the talocalcaneal axes are composed of two rotations in two cardinal planes and around two anatomical axes. To determine the axis of the talocrural articulation, a rotation around the superoinferior axis in the transverse plane must be completed (see Fig. 2a), followed by a rotation around the anteroposterior axis in the frontal plane (see Fig. 2b). For the talocalcaneal articulation, a rotation around the superoinferior axis in the transverse plane (see Fig. 2a) is followed by a rotation around the mediolateral axis in the sagittal plane (see Fig. 2c). Since the movement of a human articulation takes place around its axis, accurately modeling the second particularity of the ankle joint complex of the human foot is of utmost importance. In fact, if the axis of the articulation is incorrectly determined, the joint of the biomechanical model may allow movements that are not anatomically expected. Thus, the anatomical accuracy of the model becomes compromised, which may strongly influence the obtained results.

Fig. 2

Orientation of the talocrural and the talocalcaneal axes in the (a) transverse, (b) frontal, and (c) sagittal planes (Color figure online)

In general, the human foot and the ankle joint complex can be modeled in distinct manners by considering approaches with a varying degree of complexity, which mimic with more or less accuracy the anatomical characteristics intrinsic to the foot [2, 3, 6, 26–29]. On the one hand, some models consider only the dorsiflexion and plantarflexion movements of the human foot. This is achieved by considering that the foot articulates with the leg via the talocrural articulation only, which is formulated as an ideal revolute joint, either in two- [9, 30–32] or three-dimensions [33, 34]. On the other hand, some studies have also proposed more elaborate models considering the movements allowed by both the talocrural and the talocalcaneal articulations of the human foot. In this sense, the ankle joint complex has been modeled either as a spherical joint [35–38], as two [39–45] or three [46] separate revolute joints, as a classical universal joint [22, 47], as two equivalent one degree-of-freedom spatial parallel mechanisms [48–50] or using coupling curves [50]. Other models for the ankle joint complex, which use two separate revolute joints for the talocrural and the talocalcaneal joints, can be found in studies using musculoskeletal modeling software, such as OpenSim or Anybody [51–54]. However, modeling the ankle joint complex as a spherical or a classical universal joint precludes the possibility of considering that the talus separates the two articulations, preventing the talocrural and the talocalcaneal axes from intersecting. This difficulty can be addressed using two separate revolute joints, one for each articulation, and adding the talus bone as one segment. Nevertheless, including the talus, which has a small mass and inertia, as a rigid segment on the biomechanical model not only increases the complexity of the model but also potentially leads to numerical issues during the resolution of the equations of motion due to the calculation of high accelerations, ultimately penalizing the computational efficiency.

The motivation behind this work arises from the fact that the available studies have not yet considered sufficiently anatomically accurate and computationally efficient modeling of the ankle joint complex, which plays a key role in studying the movement of the human foot. Thus, this work aims to present a new skeletal model for the ankle joint complex of the human foot, which considers the two specific and crucial anatomical particularities described above. The proposed formulation has its foundations in a previously published work [55, 56], and it relies on the use of a modified universal joint, which considers both the talocrural and the talocalcaneal joints and is incorporated with a massless link representing the talus, ensuring the preservation of the distance between the two joint axes. The proposed approach is compared with the model presented by Anderson and Pandy [22], which is based on a classical universal joint.

The remaining of this paper is organized as follows. Section 2 briefly describes the Newton-Euler equations of motion for constrained multibody systems. The complete characterization of two different methodologies to model the ankle joint complex is presented in Sect. 3. A general methodology to impose some restrictions on the range of motion in human articulations is provided in Sect. 4. The description of the biomechanical model utilized to assess and compare the two joint models utilized to formulate the ankle joint complex is provided in Sect. 5. In the aftermath of this process, the obtained results are discussed in Sect. 6 to demonstrate the influence of each ankle joint model on the dynamic response of the human foot. Finally, this work ends with the concluding remarks in Sect. 7.

2 Multibody systems methodology

A multibody system comprises two main features, namely interconnected rigid and/or deformable bodies describing large rotational and translational displacements and joints that kinematically constrain the relative motion of the adjacent bodies. Additionally, multibody mechanical systems can also be subjected to force or driving elements [57]. Multibody dynamics involves the development of mathematical models of the systems under analysis, as well as the implementation of computational procedures to simulate, analyze, and optimize the global motion produced. The Newton-Euler approach with absolute coordinates is widely utilized to model multibody mechanical systems due to its simplicity and straightforward application to general-purpose codes [57–59]. The formulation considered in this work uses Cartesian coordinates, in which the degrees-of-freedom of the system are described by three translational ($x$, $y$, and $z$) and four orientational (Euler parameters) coordinates. Using the Baumgarte stabilization approach, the equation of motion of a general constrained multibody system can be written as

$$ \left [ \textstyle\begin{array}{c@{\quad}c} \mathbf{{M}} & \mathbf{D}^{\mathrm{T}} \\ \mathbf{D} & \mathbf{0} \end{array}\displaystyle \right ]\left \{ \textstyle\begin{array}{c} \dot{\mathbf{v}} \\ \boldsymbol{\lambda} \end{array}\displaystyle \right \} = \left \{ \textstyle\begin{array}{c} \mathbf{g} \\ \boldsymbol{\upgamma} - {2}\alpha \dot{\boldsymbol{\Phi}} - \beta ^{2}\boldsymbol{\Phi} \end{array}\displaystyle \right \} $$

(1)

in which $\mathbf{M}$ is the system mass matrix, $\mathbf{D}$ represents the Jacobian matrix of the kinematic constraint equations,$\dot{\mathbf{v}}$ denotes the vector containing the system accelerations, $\boldsymbol{\lambda}$ is the Lagrange multipliers vector associated with the reaction forces and moments on the kinematic joints and it is dependent on the formulation of the constraint equations, $\mathbf{g}$ represents the vector of applied forces and moments, and $\boldsymbol{\upgamma}$ denotes the right-hand side vector of the acceleration constraint equations. Parameters $\alpha$ and $\beta$ are the Baumgarte feedback control coefficients for velocity and position constraint violations, respectively, and are taken as positive constants. The vectors $\boldsymbol{\Phi}$ and $\dot{\boldsymbol{\Phi}}$ represent the position and velocity constraint equations, respectively [59–63]. These quantities will be explained in detail in the Sect. 3.

When the Baumgarte stabilization method is not considered in Eq. (1), the resulting equation does not explicitly include the position and velocity constraint equations. In the first steps of the simulation, the constraint violations are usually small and negligible. However, for moderate or long simulations, the violation of the original constraints formulated for the problem cannot be neglected. Thus, after several simulation steps, the results might be unacceptable [59–62, 64].

Equation (1) may be described as a combination of the differential equations of motion and the algebraic kinematic constraint equations, often referred to as a set of differential algebraic equations. The first step in the standard resolution of the equations of motion deals with the definition of the initial instant of time, $t_{0}$, and the initial conditions of the system, specifically the initial positions, $\mathbf{q}_{0}$, and velocities, $\mathbf{v}_{0}$. Subsequently, $\mathbf{M}$, $\mathbf{D}$, $\boldsymbol{\upgamma}$, $\mathbf{g}$, $\boldsymbol{\Phi}$ and $\dot{\boldsymbol{\Phi}}$ are evaluated, and Eq. (1) is solved for the constrained multibody system to obtain $\dot{\mathbf{v}}$ and $\boldsymbol{\lambda}$ at time instant $t$. Then, the auxiliary vector $\dot{\mathbf{y}}_{t}$ for time $t$, which contains the velocities and accelerations of the system, is assembled and numerically integrated for time instant $t + \Delta t$ to obtain vector $\mathbf{y}$_t + Δt. Thus, the new velocities and positions of the system are obtained. In order to proceed with the simulation, the time variable must be updated after each time step is completed. If the final simulation time has not been reached, the simulation continues, and the steps described above are repeated for the following time steps. This procedure continues until the final time of the analysis is reached [60]. The flowchart of the computational procedure for the dynamic analysis of multibody mechanical systems is presented in Fig. 3.

Fig. 3

Flowchart illustrating the computational procedure for the dynamic analysis of multibody systems based on the Baumgarte stabilization method [60]

3 Kinematic modeling of the ankle joint complex

3.1 Classical universal joint

Anderson and Pandy [22] utilized a classical universal joint to model the ankle joint complex of the human foot, as can be observed in Fig. 4a. This type of kinematic joint, also known as Hooke or Cardan joint, consists of a pair of hinge joints, the axes of which intersect at a given point $P$, and are oriented 90° relative to each other [59]. A classical universal joint allows two relative degrees-of-freedom between the bodies $i$ and $j$ that are connected by the joint. This joint transmits motion through an intermediate body, known as cross or spider, with reduced inertia and dimensions, which is not modeled per se. Figure 4b shows a general representation of a classical universal joint connecting bodies $i$ and $j$, where their centers of mass are $O_{i}$ and $O_{j}$, respectively. Body-fixed coordinate systems $\xi \eta \zeta$ are attached to the center of mass, while xyz represents the global coordinate system [59].

Fig. 4

(a) Schematic representation of the ankle joint complex of the human foot modeled as it was considered in the work of Anderson and Pandy [22]. In the anatomical reference position, the $x$, $y$, and $z$ directions correspond to the anteroposterior, mediolateral, and superoinferior directions, respectively. (b) Generic configuration of a classical universal joint connecting bodies $i$ and $j$ (Color figure online)

According to Nikravesh, the kinematic constraint equations utilized to model a classical universal joint are expressed as [59]

$$\begin{aligned} & \boldsymbol{\Phi}^{(\mathrm{s},3)} \equiv \mathbf{r}_{j}^{P} - \mathbf{r}_{i}^{P} = \mathbf{r}_{j} + \mathbf{s}_{j}^{P} - \mathbf{r}_{i} - \mathbf{s}_{i}^{P} = \mathbf{0} \end{aligned}$$

(2)

$$\begin{aligned} & {\varPhi}^{(\mathrm{n},1)} \equiv \mathbf{s}_{i}^{\mathrm{T}}\mathbf{s}_{j} = 0 \end{aligned}$$

(3)

in which $\mathbf{r}_{k}^{P}$ represents the global position vector of point $P$ located on body $k$ ($k = i$, $j$), $\mathbf{r}_{k}$ denotes the position vector of the center of mass of body $k$ described in global coordinates, and $\mathbf{s}_{k}^{P}$ is the global position vector of point $P$ located on body $k$ with respect to the body’s local coordinate system [59]. The local and global components of $\mathbf{s}_{k}^{P}$ are related to each other as

$$ \mathbf{s}_{k}^{P} = \mathbf{A}_{k}\mathbf{s}\boldsymbol{'}_{k}^{P} \left ( k = i,j \right ) $$

(4)

where $\mathbf{A}_{k}$ is the rotation matrix of body $k$, and $\mathbf{s}\boldsymbol{'}_{k}^{P}$ refers to the position of point $P$ on the body’s local coordinates. In addition, $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ represent two vectors placed on the joint axes of the classical universal joint that belong to bodies $i$ and $j$, respectively, in global coordinates. Following the Nikravesh formulation,$\boldsymbol{\Phi}^{(\mathrm{s},3)}$refers to a spherical joint constraint formulation containing three equations, and ${\varPhi}^{(\mathrm{n},1)}$ represents a normal constraint between two vectors [59].

Equation (2) establishes that the geometric center of the classical universal joint has constant coordinates with respect to any of the local coordinate systems of the connected bodies. This means that point $P_{i}$ on body $i$ must be coincident with point $P_{j}$ on body $j$, thus restricting the relative position between those bodies. In turn, Eq. (3) ensures that the two vectors $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ must remain perpendicular [59]. However, instead of considering that the axes of the ankle joint complex were perpendicular to each other, Anderson and Pandy [22] considered a specific orientation, as displayed in Fig. 2. As stated in Sect. 1, the orientation of both joint axes is composed of two rotations in two cardinal planes. In this sense, the angles shown in Fig. 2 are projected onto the corresponding plane. Therefore, after the first rotation in one plane using the angle projected on that plane, the second rotation requires calculating the corresponding angle as the real angle between the transverse plane of the body and the joint axis. In order to determine the real angles, the following conditions must be utilized [65, 66]

$$\begin{aligned} \tan \theta _{\mathrm{r}\mathrm{e}} = \tan \theta \cos \rho \end{aligned}$$

(5)

$$\begin{aligned} \tan \tau _{\mathrm{r}\mathrm{e}} = \tan \tau \cos \varepsilon \end{aligned}$$

(6)

for the talocalcaneal and talocrural articulations, respectively. The subscript ‘re’ represents the real angle. The values of $\rho$, $\theta$, $\varepsilon$ and $\tau$ are 19.8°, 35.3°, 6.0°, and 8.0°, respectively, and represent average angles retrieved from the work of Anderson and Pandy [22]. After calculating the real angles using Eqs. (5) and (6), the talocalcaneal, $\mathbf{s}_{i}$, and the talocrural, $\mathbf{s}_{j}$, axes are expressed as (see Fig. 4a)

$$\begin{aligned} & \mathbf{s}_{i} = \left \{ \textstyle\begin{array}{c@{\quad}c@{\quad}c} \cos \rho \cos \theta _{\mathrm{re}} & \sin \rho \cos \theta _{\mathrm{re}} & \sin \theta _{\mathrm{re}} \end{array}\displaystyle \right \}^{\mathrm{T}} \end{aligned}$$

(7)

$$\begin{aligned} & \mathbf{s}_{j} = \left \{ \textstyle\begin{array}{c@{\quad}c@{\quad}c} \sin \varepsilon \cos \tau _{\mathrm{re}} & \cos \varepsilon \cos \tau _{\mathrm{re}} & \sin \tau _{\mathrm{re}} \end{array}\displaystyle \right \}^{\mathrm{T}} \end{aligned}$$

(8)

In order to accommodate the fact that the axes of the ankle joint complex are non-orthogonal and have a specific orientation, Eq. (3) is modified as

$$ {\varPhi}^{(\mathrm{n},1)} \equiv \mathbf{s}_{i}^{\mathrm{T}}\mathbf{s}_{j} - \mathbf{s}_{i,0}^{\mathrm{T}}\mathbf{s}_{j,0} = 0 $$

(9)

where $\mathbf{s}$_i,0 and $\mathbf{s}$_j,0 are the coordinates of vectors $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ at the initial configuration, respectively.

With the incorporation of Eq. (9), Anderson and Pandy [22] fulfilled only one of the anatomical particularities of the ankle joint complex, neglecting that the two joint axes are non-coplanar. Thus, Eqs. (2) and (9) can be combined to establish the kinematic constraint equations of the classical universal joint utilized by Anderson and Pandy [22] as

$$ \boldsymbol{\Phi}^{(\mathrm{u},4)} \equiv \left \{ \textstyle\begin{array}{l} \mathbf{r}_{j} + \mathbf{s}_{j}^{P} - \mathbf{r}_{i} - \mathbf{s}_{i}^{P} = \mathbf{0} \hfill \\ \mathbf{s}_{i}^{\mathrm{T}}\mathbf{s}_{j} - \mathbf{s}_{i,0}^{\mathrm{T}}\mathbf{s}_{j,0} = 0 \hfill \end{array}\displaystyle \right . $$

(10)

The velocity constraint equations for a classical universal joint can be obtained by taking the first time derivative of Eq. (10), yielding

$$ \dot{\boldsymbol{\Phi}}^{(\mathrm{u},4)} \equiv \left \{ \textstyle\begin{array}{l} \dot{\mathbf{r}}_{j} + \dot{\mathbf{s}}^{P}_{j} - \dot{\mathbf{r}}_{i} - \dot{\mathbf{s}}^{P}_{i} = \mathbf{0} \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} + \mathbf{s}_{i}^{\mathrm{T}}\dot{\mathbf{s}}_{j} = 0 \hfill \end{array}\displaystyle \right . $$

(11)

in which the dot represents the derivative with respect to time.

The formulation of the classical universal joint relies exclusively on four kinematic constraints that are not explicit functions of time, as can be observed in Eq. (10). Thus, in this particular case, vector $\boldsymbol{\nu}$, which represents the right-hand side vector of the velocity constraint equations, is null.

Equation (11) is not explicitly formulated in terms of the angular velocities of bodies $i$ and $j$, which are essential for establishing the Jacobian matrix. In order to address this issue, a condition that relates the linear and angular velocities should be considered as [59]

$$ \dot{\mathbf{s}} = \tilde{\boldsymbol{\upomega}} \mathbf{s} = - \tilde{\mathbf{s}}\boldsymbol{\upomega} $$

(12)

where the symbol tilde (∼) represents the skew-symmetric matrix associated with that vector, and $\boldsymbol{\upomega}$ is the angular velocity vector in global coordinates. Thus, Eq. (11) can be rewritten as

$$ \dot{\boldsymbol{\Phi}}^{(\mathrm{u},4)} \equiv \left \{ \textstyle\begin{array}{l} \dot{\mathbf{r}}_{j} - \tilde{\mathbf{s}}_{j}^{P}\boldsymbol{\upomega}_{j} - \dot{\mathbf{r}}_{i} + \tilde{\mathbf{s}}_{i}^{P}\boldsymbol{\upomega}_{i} = \mathbf{0} \hfill \\ - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}\boldsymbol{\upomega}_{i} - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}\boldsymbol{\upomega}_{j} = 0 \hfill \end{array}\displaystyle \right . $$

(13)

or, alternatively, in the matrix form

$$ \dot{\boldsymbol{\Phi}}^{(\mathrm{u},4)} \equiv \left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} - \mathbf{I} & \tilde{\mathbf{s}}_{i}^{P} & \mathbf{I} & - \tilde{\mathbf{s}}_{j}^{P} \\ \mathbf{0} & - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i} & \mathbf{0} & - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j} \end{array}\displaystyle \right ]\left \{ \textstyle\begin{array}{c} \dot{\mathbf{r}}_{i} \\ \boldsymbol{\upomega}_{i} \\ \dot{\mathbf{r}}_{j} \\ \boldsymbol{\upomega}_{j} \end{array}\displaystyle \right \}=\mathbf{0} $$

(14)

Thus, observing Eq. (14), the contribution to the Jacobian matrix of the classical universal joint constraints can be expressed as

$$ \mathbf{D}^{(\mathrm{u},4)} = \left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} - \mathbf{I} & \tilde{\mathbf{s}}_{i}^{P} & \mathbf{I} & - \tilde{\mathbf{s}}_{j}^{P} \\ \mathbf{0} & - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i} & \mathbf{0} & - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j} \end{array}\displaystyle \right ] $$

(15)

The time derivative of Eq. (11) yields the acceleration constraint equations of the classical universal joint as

$$ \ddot{\boldsymbol{\Phi}}^{(\mathrm{u},4)} \equiv \left \{ \textstyle\begin{array}{l} \ddot{\mathbf{r}}_{j} + \ddot{\mathbf{s}}^{P}_{j} - \ddot{\mathbf{r}}_{i} - \ddot{\mathbf{s}}^{P}_{i} = \mathbf{0} \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}\ddot{\mathbf{s}}_{i} + \dot{\mathbf{s}}_{i}^{\mathrm{T}}\dot{\mathbf{s}}_{j} + \mathbf{s}_{i}^{\mathrm{T}}\ddot{\mathbf{s}}_{j} + \dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} = 0 \hfill \end{array}\displaystyle \right . $$

(16)

Equation (16) is not expressed in terms of the angular accelerations of bodies $i$ and $j$, which are required to establish the right-hand side vector of the acceleration equations. In order to address this issue, a condition relating the linear and angular accelerations must be utilized [59]

$$ \ddot{\mathbf{s}} = - \tilde{\mathbf{s}}\dot{\boldsymbol{\upomega}} + \tilde{\boldsymbol{\upomega}} \dot{\mathbf{s}} $$

(17)

where $\dot{\boldsymbol{\upomega}}$ is the angular acceleration vector in global coordinates. Thus, Eq. (16) can be rewritten as

$$ \ddot{\boldsymbol{\Phi}}^{(\mathrm{u},4)} \equiv \left \{ \textstyle\begin{array}{l} \ddot{\mathbf{r}}_{j} - \tilde{\mathbf{s}}_{j}^{P}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} + \tilde{\mathbf{s}}_{i}^{P}\dot{\boldsymbol{\upomega}}_{i} - \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P} = \mathbf{0} \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}( - \tilde{\mathbf{s}}_{i}\dot{\boldsymbol{\upomega}}_{i} + \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}) + \mathbf{s}_{i}^{\mathrm{T}}( - \tilde{\mathbf{s}}_{j}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}) + 2\dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} = 0 \hfill \end{array}\displaystyle \right . $$

(18)

By observing Eq. (18), the contribution to the right-hand side of the acceleration equations of the classical universal joint constraints can be established as

$$ \boldsymbol{\upgamma}^{(\mathrm{u},4)} = \left \{ \textstyle\begin{array}{c} \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P} - \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} \\ - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i} - 2\dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} \end{array}\displaystyle \right \} $$

(19)

3.2 Modified universal joint – proposed formulation

In this section, a modified universal joint, which is incorporated with a massless link that implicitly models the talus bone and mimics its real function, is utilized to model the ankle joint complex of the human foot, as can be observed in Fig. 5a. The proposed formulation closely follows the work developed by Malaquias et al. [55, 56] and incorporates the two anatomical particularities associated with the ankle joint complex presented in Sect. 1. A general representation of the proposed modified universal joint connecting bodies $i$ and $j$ is depicted in Fig. 5b. The centers of mass of the bodies are $O_{i}$ and $O_{j}$. Body-fixed coordinate systems $\xi \eta \zeta$ are attached to the center of mass, while xyz represents the global coordinate system [59].

Fig. 5

(a) Schematic representation of the ankle joint complex of the human foot modeled as proposed in this work. In the anatomical reference position, the $x$, $y$, and $z$ directions correspond to the anteroposterior, mediolateral, and superoinferior directions, respectively. (b) Generic configuration of a modified universal joint connecting bodies $i$ and $j$ (Color figure online)

Similar to the classical universal joint, the modified universal joint allows two degrees of freedom. Thus, four kinematic constraint equations must be considered, namely:

($i$) Establishing a constant length between points $P_{i}$ and $P_{j}$ on each of the joint axes (see Fig. 5). This constraint represents a link connecting the two bodies, which is not modeled as an actual body but instead as a massless element [57, 59]. The massless link corresponds to the talus bone of the ankle joint complex, and the constraint equation can be written as

$$ {\varPhi}_{1} \equiv \mathbf{d}^{\mathrm{T}}\mathbf{d} - l^{2} = 0 $$

(20)

where $l$ is the length of the massless link and vector $\mathbf{d}$ is defined as (see Fig. 5)

$$ \mathbf{d} = \mathbf{r}_{j}^{P} - \mathbf{r}_{i}^{P} = \mathbf{r}_{j} + \mathbf{s}_{j}^{P} - \mathbf{r}_{i} - \mathbf{s}_{i}^{P} $$

(21)

($ii$) Establishing a constant angle between the non-intersecting $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ joint axes vectors (see Fig. 5). This constraint ensures that the relative orientation of vectors $\mathbf{s}_{i}$ and $\mathbf{s}_{j}$ of the modified universal joint at any given configuration is equal to the initial one, and it is formulated as

$$ {\varPhi}_{2} \equiv \mathbf{s}_{i}^{\mathrm{T}}\mathbf{s}_{j} - \mathbf{s}_{i,0}^{\mathrm{T}}\mathbf{s}_{j,0} = 0 $$

(22)

(iii) Establishing a constant angle between vector $\mathbf{d}$ and the $\mathbf{s}_{i}$ joint axis vector (see Fig. 5). This constraint is added to the formulation of the modified universal joint to maintain the relative orientation of the joint axis $\mathbf{s}_{i}$ and vector $\mathbf{d}$ constant in all configurations. The corresponding kinematic constraint equation can be expressed as

$$ {\varPhi}_{3} \equiv \mathbf{s}_{i}^{\mathrm{T}}\mathbf{d} - \mathbf{s}_{i,0}^{\mathrm{T}}\mathbf{d}_{0} = 0 $$

(23)

where $\mathbf{d}_{0}$ represents the initial coordinates of vector $\mathbf{d}$.

($iv$) Establishing a constant angle between vector $\mathbf{d}$ and the $\mathbf{s}_{j}$ joint axis vector (see Fig. 5). The constraint is expressed as

$$ {\varPhi}_{4} \equiv \mathbf{s}_{j}^{\mathrm{T}}\mathbf{d} - \mathbf{s}_{j,0}^{\mathrm{T}}\mathbf{d}_{0} = 0 $$

(24)

In summary, the set of kinematic constraint equations utilized to model the modified universal joint proposed in this work can be condensed as

$$ \boldsymbol{\Phi}^{(\mathrm{u}_{\mathrm{m}},4)} \equiv \left \{ \textstyle\begin{array}{l} \mathbf{d}^{\mathrm{T}}\mathbf{d} - l^{2} = 0 \hfill \\ \mathbf{s}_{i}^{\mathrm{T}}\mathbf{s}_{j} - \mathbf{s}_{i,0}^{\mathrm{T}}\mathbf{s}_{j,0} = 0 \hfill \\ \mathbf{s}_{i}^{\mathrm{T}}\mathbf{d} - \mathbf{s}_{i,0}^{\mathrm{T}}\mathbf{d}_{0} = 0 \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}\mathbf{d} - \mathbf{s}_{j,0}^{\mathrm{T}}\mathbf{d}_{0} = 0 \hfill \end{array}\displaystyle \right . $$

(25)

The velocity constraint equations for the proposed modified universal joint are obtained by taking the first time derivative of Eq. (25) yielding

$$ \dot{\boldsymbol{\Phi}}^{(\mathrm{u}_{\mathrm{m}},4)} \equiv \left \{ \textstyle\begin{array}{l} 2\mathbf{d}^{\mathrm{T}}(\dot{\mathbf{r}}_{j} + \dot{\mathbf{s}}_{j}^{P} - \dot{\mathbf{r}}_{i} - \dot{\mathbf{s}}_{i}^{P}) = 0 \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} + \mathbf{s}_{i}^{\mathrm{T}}\dot{\mathbf{s}}_{j} = 0 \hfill \\ \mathbf{d}^{\mathrm{T}}\dot{\mathbf{s}}_{i} + \mathbf{s}_{i}^{\mathrm{T}}(\dot{\mathbf{r}}_{j} + \dot{\mathbf{s}}_{j}^{P} - \dot{\mathbf{r}}_{i} - \dot{\mathbf{s}}_{i}^{P}) = 0 \hfill \\ \mathbf{d}^{\mathrm{T}}\dot{\mathbf{s}}_{j} + \mathbf{s}_{j}^{\mathrm{T}}(\dot{\mathbf{r}}_{j} + \dot{\mathbf{s}}_{j}^{P} - \dot{\mathbf{r}}_{i} - \dot{\mathbf{s}}_{i}^{P}) = 0 \hfill \end{array}\displaystyle \right . $$

(26)

Since the kinematic constraints of Eq. (25) are not explicit functions of time, vector $\boldsymbol{\nu}$ is null.

In a similar manner to the classical universal joint, considering the condition expressed in Eq. (12), then Eq. (26) can be rewritten as

$$ \dot{\boldsymbol{\Phi}}^{(\mathrm{u}_{\mathrm{m}},4)} \equiv \left \{ \textstyle\begin{array}{l} 2\mathbf{d}^{\mathrm{T}}(\dot{\mathbf{r}}_{j} - \tilde{\mathbf{s}}_{j}^{P}\boldsymbol{\upomega}_{j} - \dot{\mathbf{r}}_{i} + \tilde{\mathbf{s}}_{i}^{P}\boldsymbol{\upomega}_{i}) = 0 \hfill \\ - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}\boldsymbol{\upomega}_{i} - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}\boldsymbol{\upomega}_{j} = 0 \hfill \\ - \mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}\boldsymbol{\upomega}_{i} + \mathbf{s}_{i}^{\mathrm{T}}\dot{\mathbf{r}}_{j} - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P}\boldsymbol{\upomega}_{j} - \mathbf{s}_{i}^{\mathrm{T}}\dot{\mathbf{r}}_{i} + \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P}\boldsymbol{\upomega}_{i} = 0 \hfill \\ - \mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}\boldsymbol{\upomega}_{j} + \mathbf{s}_{j}^{\mathrm{T}}\dot{\mathbf{r}}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P}\boldsymbol{\upomega}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\dot{\mathbf{r}}_{i} + \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P}\boldsymbol{\upomega}_{i} = 0 \hfill \end{array}\displaystyle \right . $$

(27)

or, alternatively, in the matrix form

$$ \dot{\boldsymbol{\Phi}}^{(\mathrm{u}_{\mathrm{m}},4)} \equiv \left [ \textstyle\begin{array}{c} - 2\mathbf{d}^{\mathrm{T}} \\ \mathbf{0} \\ - \mathbf{s}_{i}^{\mathrm{T}} \\ - \mathbf{s}_{j}^{\mathrm{T}} \end{array}\displaystyle \textstyle\begin{array}{c} \\ \\ \\ \end{array}\displaystyle \textstyle\begin{array}{c} 2\mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P} \\ - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i} \\ - \mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{i} + \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P} \\ \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P} \end{array}\displaystyle \textstyle\begin{array}{c} \\ \\ \\ \end{array}\displaystyle \textstyle\begin{array}{c} 2\mathbf{d}^{\mathrm{T}} \\ \mathbf{0} \\ \mathbf{s}_{i}^{\mathrm{T}} \\ \mathbf{s}_{j}^{\mathrm{T}} \end{array}\displaystyle \textstyle\begin{array}{c} \\ \\ \\ \end{array}\displaystyle \textstyle\begin{array}{c} - 2\mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P} \\ - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j} \\ - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P} \\ - \mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P} \end{array}\displaystyle \right ]\left \{ \textstyle\begin{array}{c} \dot{\mathbf{r}}_{i} \\ \boldsymbol{\upomega}_{i} \\ \dot{\mathbf{r}}_{j} \\ \boldsymbol{\upomega}_{j} \end{array}\displaystyle \right \} =\mathbf{0} $$

(28)

Thus, observing Eq. (28), the contribution of the modified universal joint constraints to the Jacobian matrix can be established as

$$ \mathbf{D}_{(4 \times 12)}^{{(\mathrm{u}_{\mathrm{m}},4)}} = \left [ \textstyle\begin{array}{c} - 2\mathbf{d}^{\mathrm{T}} \\ \mathbf{0} \\ - \mathbf{s}_{i}^{\mathrm{T}} \\ - \mathbf{s}_{j}^{\mathrm{T}} \end{array}\displaystyle \textstyle\begin{array}{c} \\ \\ \\ \end{array}\displaystyle \textstyle\begin{array}{c} 2\mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P} \\ - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i} \\ - \mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{i} + \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P} \\ \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{i}^{P} \end{array}\displaystyle \textstyle\begin{array}{c} \\ \\ \\ \end{array}\displaystyle \textstyle\begin{array}{c} 2\mathbf{d}^{\mathrm{T}} \\ \mathbf{0} \\ \mathbf{s}_{i}^{\mathrm{T}} \\ \mathbf{s}_{j}^{\mathrm{T}} \end{array}\displaystyle \textstyle\begin{array}{c} \\ \\ \\ \end{array}\displaystyle \textstyle\begin{array}{c} - 2\mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P} \\ - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j} \\ - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P} \\ - \mathbf{d}^{\mathrm{T}}\tilde{\mathbf{s}}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\mathbf{s}}_{j}^{P} \end{array}\displaystyle \right ] $$

(29)

Finally, the time derivative of Eq. (26) yields the acceleration constraint equations of the modified universal joint as

$$ \ddot{\boldsymbol{\Phi}}^{(\mathrm{u}_{\mathrm{m}},4)} \equiv \left \{ \textstyle\begin{array}{l} 2\mathbf{d}^{\mathrm{T}}(\ddot{\mathbf{r}}_{j} + \ddot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} - \ddot{\mathbf{s}}_{i}^{P}) + 2\dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{d}} = 0 \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}\ddot{\mathbf{s}}_{i} + \dot{\mathbf{s}}_{i}^{\mathrm{T}}\dot{\mathbf{s}}_{j} + \mathbf{s}_{i}^{\mathrm{T}}\ddot{\mathbf{s}}_{j} + \dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} = 0 \hfill \\ \mathbf{d}^{\mathrm{T}}\ddot{\mathbf{s}}_{i} + \dot{\mathbf{s}}_{i}^{\mathrm{T}}\dot{\mathbf{d}} + \mathbf{s}_{i}^{\mathrm{T}}(\ddot{\mathbf{r}}_{j} + \ddot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} - \ddot{\mathbf{s}}_{i}^{P}) + \dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{s}}_{i} = 0 \hfill \\ \mathbf{d}^{\mathrm{T}}\ddot{\mathbf{s}}_{j} + \dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{d}} + \mathbf{s}_{j}^{\mathrm{T}}(\ddot{\mathbf{r}}_{j} + \ddot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} - \ddot{\mathbf{s}}_{i}^{P}) + \dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{s}}_{j} = 0 \hfill \end{array}\displaystyle \right . $$

(30)

Taking advantage of Eq. (17), then Eq. (30) can be rewritten as

$$ \ddot{\boldsymbol{\Phi}}^{(\mathrm{u}_{\mathrm{m}},4)} \equiv \left \{ \textstyle\begin{array}{l} 2\mathbf{d}^{\mathrm{T}}(\ddot{\mathbf{r}}_{j} - \tilde{\mathbf{s}}_{j}^{P}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} + \tilde{\mathbf{s}}_{i}^{P}\dot{\boldsymbol{\upomega}}_{i} - \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P}) + 2\dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{d}} = 0 \hfill \\ \mathbf{s}_{j}^{\mathrm{T}}( - \tilde{\mathbf{s}}_{i}\dot{\boldsymbol{\upomega}}_{i} + \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}) + \dot{\mathbf{s}}_{i}^{\mathrm{T}}\dot{\mathbf{s}}_{j} + \mathbf{s}_{i}^{\mathrm{T}}( - \tilde{\mathbf{s}}_{j}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}) + \dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} = 0 \hfill \\ \mathbf{d}^{\mathrm{T}}( - \tilde{\mathbf{s}}_{i}\dot{\boldsymbol{\upomega}}_{i} + \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}) + \dot{\mathbf{s}}_{i}^{\mathrm{T}}\dot{\mathbf{d}} + \mathbf{s}_{i}^{\mathrm{T}}(\ddot{\mathbf{r}}_{j} - \tilde{\mathbf{s}}_{j}^{P}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} + \tilde{\mathbf{s}}_{i}^{P}\dot{\boldsymbol{\upomega}}_{i} - \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P}) + \dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{s}}_{i} = 0 \hfill \\ \mathbf{d}^{\mathrm{T}}( - \tilde{\mathbf{s}}_{j}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}) + \dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{d}} + \mathbf{s}_{j}^{\mathrm{T}}(\ddot{\mathbf{r}}_{j} - \tilde{\mathbf{s}}_{j}^{P}\dot{\boldsymbol{\upomega}}_{j} + \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} - \ddot{\mathbf{r}}_{i} + \tilde{\mathbf{s}}_{i}^{P}\dot{\boldsymbol{\upomega}}_{i} - \tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P}) + \dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{s}}_{j} = 0 \hfill \end{array}\displaystyle \right . $$

(31)

Thus, the contribution of the modified universal joint constraints to the right-hand side of the acceleration equations yields

$$ \boldsymbol{\upgamma}^{{(\mathrm{u}_{\mathrm{m}},4)}} = \left \{ \textstyle\begin{array}{c} 2\mathbf{d}^{\mathrm{T}}(\tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P} - \tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P}) - 2\dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{d}} \\ - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i} - 2\dot{\mathbf{s}}_{j}^{\mathrm{T}}\dot{\mathbf{s}}_{i} \\ - \mathbf{d}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i} - 2\dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{s}}_{i} - \mathbf{s}_{i}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} + \mathbf{s}_{i}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P} \\ - \mathbf{d}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j} - 2\dot{\mathbf{d}}^{\mathrm{T}}\dot{\mathbf{s}}_{j} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{j}\dot{\mathbf{s}}_{j}^{P} - \mathbf{s}_{j}^{\mathrm{T}}\tilde{\boldsymbol{\upomega}}_{i}\dot{\mathbf{s}}_{i}^{P} \end{array}\displaystyle \right \} $$

(32)

In comparison with the ankle joint complex models presented in Sect. 1, the use of the massless link approach has several advantages, namely ($i$) it is guaranteed that the physiological distance between the talocrural and the talocalcaneal articulations is preserved; ($ii$) the specific anatomical orientations of the joint axes are taken into account, (iii) the impossibility of measuring the talus motion using non-invasive techniques is avoided, ($iv$) the overall complexity of the biomechanical model is not increased as the number of coordinates and constraints is kept unchanged, and ($v$) the addition of a segment with very small mass and inertia in the mass matrix is avoided. This last issue is of paramount importance in the measure that a segment with small mass and inertia would result in the calculation of significantly high accelerations during the resolution of the equations of motion defined in Eq. (1). The calculation of such high values for the accelerations would ultimately decrease the integration time of the system and lower its computational efficiency. Concerning the advantage (iii), since the talus is an internal bone that is hidden by the surrounding tibia, fibula, and calcaneus bones, its motion is impossible to measure when acquiring experimental data of the human movement with non-invasive techniques, such as the use of skin markers [45, 67]. Introducing the massless link in the proposed formulation overcomes this difficulty, eliminating the need for experimental measurement of talus motion.

4 Methodology to restrict the range of motion in human articulations

This section deals with a methodology to restrict the range of motion in human articulations. The formulation presented here closely follows the work developed by Silva et al. [68], and it can be applied to any general rotational joint.

The amplitude of motion allowed in any anatomical articulation is usually known as the range of motion, and it is strongly influenced by the bony structures adjacent to the articulation, as well as by the physiological properties of the surrounding muscles and ligaments. Generally, two distinct types of range of motion can be considered, namely normal and limited. The normal range of motion is associated with healthy articulations, while the limited range of motion is related to pathological conditions or injuries. It is clear that the normal range of motion allows the articulation to adjust easily to perturbations imposed by the human body, reducing the risk of injury. In turn, a limited range of motion impairs function, mobility, and daily life activities, eventually leading to pain [69–72]. It is also important to distinguish between active and passive ranges of motion. Both the active and the passive ranges of motion allow local movement of the articulations. The active range of motion occurs when the muscles contract and relax, while the passive range of motion occurs when external loads act on the articulation. The passive range of motion is usually the maximum range of motion of the articulation, and it is limited by pain and by the resistance due to the tissues existing around the articulation [23, 73, 74].

From the multibody systems formulation point of view, the kinematic joint models presented in the previous section do not impose any restrictions on the amplitude of motion permitted between the adjacent bodies connected by them. Therefore, it is necessary to incorporate additional restrictions to prevent the joints from performing anatomically and physiologically unacceptable movements. In what follows, an approach to restrict the range of motion in human articulations is presented under the umbrella of multibody systems formulation.

Joint resistance moments, which mimic the passive behavior of the tissues surrounding human articulations [75–77], are applied to prevent unrealistic configurations of the articulations during regular and safe daily activities. The joint resistance moments associated with a generic rotational joint $h$ can be established as [68]

$$ \mathbf{m}_{h}^{\mathrm{r}} = \mathbf{m}_{h}^{\mathrm{d}} + \mathbf{m}_{h}^{\mathrm{mr}} $$

(33)

where $\mathbf{m}_{h}^{\mathrm{d}}$ is the dissipative moment vector and $\mathbf{m}_{h}^{\mathrm{mr}}$ represents the motion-restricting moment vector. The subscript $h$ denotes de the joint.

In order to better understand the presented methodology, it is important to clarify the concepts of reference and moving bodies. A schematic representation of bodies $i$ and $j$, which are, respectively, the reference and moving bodies, connected by a generic joint $h$ is presented in Fig. 6. While the moving body is the one whose range of motion is intended to be restricted, the reference body enables the definition of the limits inside which the moving body can freely move without exceeding the joint’s range of motion. It must be noticed that the joint resistance moments calculated using Eq. (33) are applied on the moving body, and their symmetrical counterparts are applied on the reference body. In the multibody systems methodology, these moments are treated as external moments introduced in vector $\mathbf{g}$ of Eq. (1).

Fig. 6

Schematic representation of a system composed of bodies $i$ and $j$, which represent, respectively, the reference and moving bodies connected by a general rotational joint $h$. The circumduction cone, the joint’s local reference $\xi _{h}\eta _{h}\zeta _{h}$ and the unit vector $\mathbf{u}$_r,h are represented (Color figure online)

The joint dissipative moment can be calculated using a viscous torsional damper as

$$ \mathbf{m}_{h}^{\mathrm{d}} = - j_{h}\boldsymbol{\upomega}_{h} $$

(34)

in which $j_{h}$ is the rotational damping coefficient of the human articulation, and $\boldsymbol{\upomega}_{h}$ represents the joint $h$ angular velocity vector. The joint dissipative moment given by Eq. (34) represents the energy lost due to the viscoelasticity of the tissues existing around human articulations.

The joint motion-restricting moment acts to restrict the range of motion of the human articulation, and it is typically modeled as a nonlinear elastic element. The joint motion-restricting moment is null for normal joint motions, that is, within the allowable range of motion for that specific joint, and increases rapidly from zero until a maximum value within a given angular range for unacceptable joint configurations. It should be noticed that the transition from null to the maximum moment magnitude occurs in a very short angular range. In the presented methodology, the acceptable configuration of a generic human articulation is established by the circumduction cone. This cone can be defined as a three-dimensional surface inside which any configuration of the articulation is possible and physiologically acceptable. Thus, the range of motion of the human articulations can be associated with the working space of the circumduction cone [68, 78], as represented in Fig. 6. The circumduction cone can be established for any joint with rotational degrees-of-freedom, independently of the number of degrees-of-freedom. Thus, for revolute joints, there are sections of the cone that are not reached by the moving body due to the kinematic structure of this type of joint. Furthermore, due to the fact that spherical and classical universal joints allow simultaneous rotation around more than one joint axis, their range of motion can be restricted by defining only one circumduction cone. It should be noted that the translational and the internal rotational degrees-of-freedom are not restricted by this methodology. Thus, only the bending moment is considered, and the torsional moment is neglected. The internal rotation has no influence on the bending moment.

In order to define the circumduction cone, the local reference frame of the reference body must be considered, which can be established in an arbitrary manner. The circumduction cone is described by defining the joint’s local reference frame, which is rigidly attached to the reference body. Thus, the axes of joint’s local reference frame are defined using the reference body’s local reference frame. It is important to note that the joint’s and reference body’s reference frames are not necessarily coincident. In the cone illustrated in Fig. 6, the unit vectors defining the axes $\xi$ and $\eta$ of the joint’s local reference frame are established as

$$ \mathbf{u}\boldsymbol{'}_{\xi ,h} = \left \{ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 0 & 0 & - 1 \end{array}\displaystyle \right \}^{\mathrm{T}}\quad \text{and}\quad \mathbf{u}\boldsymbol{'}_{\eta ,h} = \left \{ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 0 & 1 & 0 \end{array}\displaystyle \right \}^{\mathrm{T}} $$

The $\mathbf{u}_{\zeta,h}$ unit vector of the joint’s local reference frame is obtained as the cross product of the $\mathbf{u}_{\xi,h}$ and $\mathbf{u}_{\eta,h}$ unit vectors as

$$ \mathbf{u}\boldsymbol{'}_{\zeta ,h} = \tilde{\mathbf{u}}\boldsymbol{'}_{\xi ,h}\mathbf{u}\boldsymbol{'}_{\eta ,h} $$

(35)

Subsequently to the definition of the joint’s local reference frame, a unit vector, $\mathbf{u}$_r,h, must be established using the moving body’s local reference frame. This unit vector allows to know, at every configuration, the orientation of the moving body in the joint’s local reference frame, except for the internal rotation. Analyzing the cone illustrated in Fig. 6, vector $\mathbf{u}$_r,h is defined as

$$ \mathbf{u}\boldsymbol{'}_{\mathrm{r},h} = \left \{ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 1 & 0 & 0 \end{array}\displaystyle \right \}^{\mathrm{T}} $$

In order to determine whether the actual position of the moving body is inside (acceptable configuration) or outside (unacceptable configuration) the circumduction cone, the longitude, $\psi _{h}$, and latitude, $\sigma _{h}$, of the unit vector $\mathbf{u}$_r,h expressed in the joint’s local reference frame must be determined. The longitude is established by the joint’s local axis $\mathbf{u}_{\xi,\textit{h}}$ and the projection of vector $\mathbf{u}$_r,h in the plane $\xi _{h}\eta _{h}$ of the joint’s local reference frame (see Fig. 7a). Therefore, $\psi _{h}$ is calculated as

$$ \psi _{h} = \arctan \left ( \frac{ u''_{\mathrm{r},h,\eta}}{ u''_{\mathrm{r},h,\xi}} \right ) $$

(36)

where the symbol (^′′) represents the components of the unit vector $\mathbf{u}$_r,h in the joint’s local reference frame. It is important to note that, in Eq. (36), the atan2 function, which denotes the four-quadrant inverse tangent, must be utilized. The atan2 function returns values in the closed interval [$- \uppi$, $\uppi$], and, therefore, whenever $\psi _{h}$ calculated using Eq. (36) is negative, a value of 2$\uppi$ must be added to $\psi _{h}$ in order to obtain values in the closed interval [0, 2$\uppi$].

Fig. 7

Schematic representation of the (a) longitude and (b) latitude utilized in the calculation of the joint motion-restricting moment. The longitude is defined in the plane $\xi _{h}\eta _{h}$ of the joint’s local reference frame (top view of the circumduction cone) (Color figure online)

In addition, the latitude, $\sigma _{h}$, is established by the joint’s local unit vector $\mathbf{u}$_ζ,h and the unit vector $\mathbf{u}$_r,h expressed in the joint’s local reference frame (see Fig. 7b), and it is determined by

$$\begin{aligned} \sigma_{h} = \arccos \left( \mathbf{u}\boldsymbol{''}_{\mathrm{r},h} {}^{\mathrm{T}}\mathbf{u} \boldsymbol{''}_{\zeta,h} \right) \end{aligned}$$

(37)

It is important to note that vector $\mathbf{u}$_r,h can be defined in the most convenient manner depending on the characteristics of the multibody system under analysis. However, it should not be defined in such a way that its initial configuration is outside the circumduction cone, thus being in an unacceptable configuration. Although it is not mandatory, in this example, the $\mathbf{u}$_r,h unit vector and the joint’s local reference frame are defined such that the unit vectors $\mathbf{u}$_r,h, and $\mathbf{u}_{\zeta,\textit{h}}$ are coincident at the initial configuration (see Fig. 6). Thus, in this situation, the initial latitude is zero, which avoids an unacceptable configuration of the moving body and prevents the joint’s range of motion from being exceeded at the initial configuration.

The circumduction cone of a certain joint is defined by specifying the maximum allowable latitude, $\sigma$_h,max, for certain values of longitude, $\psi _{h}$. The longitude and latitude are directly associated with the joint’s range of motion and must be defined in the closed interval [0, 2$\uppi$] and [0, $\uppi$] radians, respectively. It is important to notice that the definition of the axes of the joint’s local reference frame influences the definition of the latitude and longitude values. Then, an interpolation function determines the maximum allowable latitude corresponding to any longitude value. If, for the current longitude, $\psi _{h}$, the current latitude, $\sigma _{h}$, exceeds the corresponding maximum allowable latitude, $\sigma$_h,max, then an unacceptable joint position is verified, and the joint motion-restricting moment must be applied. This moment corresponds to a third-degree polynomial function with the behavior depicted in Fig. 8, and it can be expressed as

$$ \mathbf{m}_{h}^{\mathrm{mr}} = \left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} \mathbf{0} \hfill & \text{if } \sigma _{h} \le \sigma _{h,\max } \hfill \\ m_{\mathrm{p},h}\left [ 3\left ( \frac{\kappa}{\Delta \sigma _{h}} \right )^{2} - 2\left ( \frac{\kappa}{\Delta \sigma _{h}} \right )^{3} \right ]\tilde{\mathbf{u}}_{\mathrm{r}, h}\mathbf{n} \hfill & \text{if } \sigma _{h,\max } < \sigma _{h}\ \text{and}\ \kappa \le \Delta \sigma _{h} \hfill \\ m_{\mathrm{p},h}\tilde{\mathbf{u}}_{\mathrm{r}, h}\mathbf{n} \hfill & \text{if } \kappa > \Delta \sigma _{h} \hfill \end{array}\displaystyle \right . $$

(38)

where $m$_p,h is the magnitude of the maximum penalty moment applied to restrict the joint’s motion, $\Delta \sigma _{h}$ denotes the value corresponding to the difference between the passive and active ranges of motion of joint $h$, and it allows to adjust the joint’s stiffness in the model presented in Eq. (38), $\kappa$ denotes the relative angular motion between the moving body and the limits of the circumduction cone, and it allows to quantify the angular displacement that is exceeded by the moving body with respect to the range of motion of the joint, and $\mathbf{n}$ represents the direction normal to the surface of the circumduction cone for the corresponding value of $\sigma$_h,max. It is important to note that, depending on the considered human articulation, the value of $\Delta \sigma _{h}$ varies due to the fact that both the active and passive ranges of motion vary from human articulation to human articulation. If the value of $\Delta \sigma _{h}$ is very small, the joint’s stiffness is very high, and vice-versa.

Fig. 8

Nonlinear behavior of the joint motion-restricting moment [68]

The expression presented in Eq. (38) can easily be substituted by others existing in the literature, such as the ones given in the works of Yamaguchi [77] and Nasr et al. [79], while keeping unchanged the remaining methodology described in the previous sections.

5 Application example: a biomechanical model of the right human foot and leg

A three-dimensional biomechanical multibody model implemented using an in-house code developed in Matlab was considered to compare the two ankle joint models presented in Sect. 3. The biomechanical model was utilized for forward dynamic analysis, and it is composed of three rigid bodies, namely the toes, the main foot, and the leg. The leg represents the tibia and fibula, the main foot is constituted by the tarsus and metatarsus, and the toes encompass the phalanges. The bodies of this biomechanical model are kinematically connected to each other by means of one revolute joint, connecting the toes to the main foot and representing the metatarsophalangeal articulations, and one universal joint, either classical or modified, connecting the main foot to the leg and representing the ankle joint complex. A fixed joint is applied at the center of mass of the leg segment, restraining the translation and rotation of this body. In addition, a constant angle kinematic constraint, similar to Eq. (22), is utilized on the metatarsophalangeal joint to lock the joint’s angle. This procedure was performed to focus the analysis of the results exclusively on the differences between the joint models considered for the ankle joint complex. The developed biomechanical model has three degrees of freedom, but, since the metatarsophalangeal joint was subjected to a constant angle kinematic constraint, in reality, the model has only two functional degrees of freedom, which arise from the universal joint. The numbers of each body and their corresponding local coordinate systems are shown in Fig. 9.

Fig. 9

Schematic representation of the biomechanical model of the right human foot and leg with the ankle joint complex modeled as a (a) classical and (b) modified universal joint. The body segments are represented in different colors (Color figure online)

The generic configuration of the biomechanical model is displayed in Fig. 9 and the the corresponding initial conditions are listed in Table 1. The location of the center of mass of each segment was retrieved from the work of Anderson and Pandy [22]. It was assumed that the locations of all joints and the center of mass of all segments were aligned with the midline of the foot, therefore resulting in a null $y$-coordinate. The $x$, $y$, and $z$ directions correspond to the anteroposterior, mediolateral, and superoinferior directions in the three cardinal planes, respectively. For both biomechanical models presented in Fig. 9, the location of the talocrural joint is the same, as well as the orientation of its axis. In addition, the talocalcaneal joint axis has the same orientation in both models. The main difference between the models presented in Fig. 9 is the location of the talocalcaneal joint.

Table 1 Initial conditions for the biomechanical multibody model

The dimensions and inertial properties of each body of the biomechanical model are listed in Table 2. It is important to note that the local reference frames of the bodies are located at their center of mass (see Fig. 9) and are aligned with the principal axes of inertia. The model corresponds to a 71 kg and 1.77 m male subject [22]. In addition, the total foot length was considered to be 0.27 m [22] and the length of the talus, which represents the length $l$ of the massless link, was taken as 0.0417 m [80].

Table 2 Dimensions and inertial properties of the bodies of the biomechanical multibody model

The simulation parameters utilized in all dynamic simulations and in the numerical methods required to solve the system dynamics are displayed in Table 3 [81, 82].

Table 3 Parameters used in the dynamic simulation of the biomechanical multibody model

Talocrural maximum range of motion for dorsiflexion and plantarflexion was considered to be 20° and 35°, respectively. Talocalcaneal maximum inversion and eversion were defined as 35° and 25°, respectively. These values refer to active ranges of motion, which are based on in vivo data [83]. The values for $\Delta \sigma _{h}$ and $m$_p,h were retrieved from the work of Silva et al. [68, 84] for the knee joint as 11.5° and 226 N$\boldsymbol{\cdot}$m, respectively. The damping coefficient of the ankle joint complex, $j_{h}$, was obtained by multiplying the total mass of the subject by the normalized value 0.008 N$\boldsymbol{\cdot}$m$\boldsymbol{\cdot}$s/kg/rad reported in [85, 86]. A value of 0.5680 N$\boldsymbol{\cdot}$m$\boldsymbol{\cdot}$s was utilized.

It should be noted that, in this work, the three-dimensional biomechanical multibody model is considered for forward dynamics analysis. However, when experimental data of the human movement acquired at the laboratory is utilized, inverse dynamics analysis must be considered. In this situation, the masses, lengths, moments of inertia, Euler parameters, and location and orientation of the local reference system of each body must be adjusted based on the experimentally acquired data.

6 Results and discussion

In the first dynamic analysis, the biomechanical model is released from the initial position with null velocities and under the action of the gravitational force only, which is assumed to act in the negative $z$-direction. The range of motion of the talocrural and the talocalcaneal joints is not restricted. The trajectory of the center of mass of the main foot segment in the sagittal ($x$-$z$ plane), transverse ($x$-$y$ plane), and frontal ($y$-$z$ plane) planes is represented in the plots of Fig. 10. From the analysis of Fig. 10a, it can be observed that there are no substantial differences in the type and amplitude of motion performed by the two ankle joint models in the sagittal plane. This outcome is expected because the talocrural joint, which is the joint most responsible for the movement in the sagittal plane, has the same location, and its joint axis has the same orientation on both the classical and modified universal joint models. Although the talocalcaneal joint has the same orientation in both joint models, its location differs, as can be observed in Figs. 4, 5, and 9. Thus, in the transverse (Fig. 10b) and frontal (Fig. 10c) planes, the observed differences in the trajectory of the center of mass of the main foot segment are considerable, mainly in terms of amplitude of motion and direction of movement. It should be noted that the amplitude of motion in the mediolateral direction ($y$-direction) allowed by the proposed approach is greater than the one allowed by the Anderson and Pandy’s model [22], which is observed by the fact that the blue plot is wider than the red plot in Fig. 10b and Fig. 10c. In addition, it can be noticed that the foot tends to move towards the lateral direction (negative $y$) in the Anderson and Pandy’s model [22], as opposed to the proposed approach, in which the foot tends to move towards the medial direction (positive $y$). In order to better understand these differences in the mediolateral direction, it is important to note that the main distinction between the Anderson and Pandy’s model [22] and the approach proposed in this work is the location of the talocalcaneal joint. In the Anderson and Pandy’s model [22], the talocrural and the talocalcaneal joints coincide at point $P$, corresponding to the anatomical location of the talocrural articulation. In the model proposed in this work, the two joint axes are non-coplanar due to the consideration of the massless link, which represents the talus bone. Thus, the talocalcaneal joint is closer to the center of mass of the main foot segment than in the Anderson and Pandy’s model [22]. The fact that the center of mass of the main foot is displaced from the location of the talocalcaneal joint in both joint models generates a moment of force. The magnitude of this moment is smaller in the proposed approach than in the Anderson and Pandy’s model [22] because the lever arm is also smaller. Consequently, the amplitude of motion allowed by the proposed approach in the mediolateral direction is greater than that observed for the Anderson and Pandy’s model [22]. The incorporation of the massless link explains the differences between the two ankle joint models. It is also important to note that the range of motion for talocrural plantarflexion is violated in both joint models, which is expected because the methodology to restrict the range of motion, presented in Sect. 4, is not considered in this analysis. To sum up, a comparison of the amplitude of motion and type of movement performed by the ankle two joint models is shown in the snapshots of Fig. 11.

Fig. 10

Trajectory of the center of mass of the main foot segment allowed by the two joint models considered in this work to model the ankle joint complex of the human foot in the (a) sagittal, (b) transverse and (c) frontal planes. The green marker represents the initial position of the center of mass of the main foot segment. The range of motion of the talocrural and the talocalcaneal joints is not restricted. The $x$, $y$, and $z$ directions correspond to the anteroposterior, mediolateral, and superoinferior directions in the three cardinal planes, respectively (Color figure online)

Fig. 11

Snapshots of the movement of the biomechanical model for the first second of simulation using (a) the Anderson and Pandy [22] classical universal joint (See Online Resource 1 – 11044_2023_9955_MOESM1_ESM) and (b) the proposed modified universal joint (See Online Resource 2 – 11044_2023_9955_MOESM2_ESM). The model is subjected to the gravitational action only. The range of motion of the talocrural and the talocalcaneal joints is not restricted (Color figure online)

The time evolution of the velocity and acceleration of the center of mass of the main foot segment is depicted in Fig. 12. The most important differences between both ankle joint models are observed in the mediolateral direction (Fig. 12c and Fig. 12d). The plots differ in magnitude, in which the proposed approach allows higher velocity and acceleration of the center of mass of the main foot, and they are inverted, meaning that when a positive velocity/acceleration value is observed in one joint model, the other model presents a negative value, and vice-versa. This phenomenon is associated with the talocalcaneal joint location, which is influenced by the consideration or not of the massless link. The results are quite similar concerning the anteroposterior (Fig. 12a and Fig. 12b) and superoinferior (Fig. 12e and Fig. 12f) directions. Small differences can be visible in the peaks for the velocity and acceleration, with the Anderson and Pandy’s model [22] exhibiting higher values.

Fig. 12

Influence of the joint model considered for the ankle joint complex of the human foot on the response of the biomechanical model. Velocity (a), (c), (e), and acceleration (b), (d), (f) of the center of mass of the main foot segment (Color figure online)

Figure 13 shows the phase portraits of the two ankle joint models, which include the $x$-, $y$- and $z$-components of the position of the center of mass of the main foot versus its velocity and the velocity versus the acceleration. The proposed joint model exhibits a less smooth behavior when compared with the Anderson and Pandy’s model [22], which is seen from the fact that the blue plot tends to be more scattered than the red plot. In the $x$- and $z$-directions, overall, the results show no substantial differences. However, the most substantial differences are observed in the $y$-direction, which agrees with the previous results (Fig. 10 and Fig. 12). This phenomenon is again associated with the difference in the location of the talocalcaneal joint.

Fig. 13

Phase portraits of the Anderson and Pandy’s and proposed joint models: (a) $x$-position vs. $x$-velocity, (b) $x$-velocity vs. $x$-acceleration, (c) $y$-position vs. $y$-velocity, (d) $y$-velocity vs. $y$-acceleration, (e) $z$-position vs. $z$-velocity and (f) $z$-velocity vs. $z$-acceleration (Color figure online)

A second analysis is performed, in which the biomechanical model is also released from the initial position with null velocities and under the action of the gravitational force. In addition, joint resistance moments are applied to restrict the range of motion of the talocrural and talocalcaneal joints. This analysis aims to compare the type of movement and the amplitude of motion performed by the two ankle models using a more realistic application example. The trajectory of the center of mass of the main foot is subjected to analysis, and it can be observed in the three cardinal planes in Fig. 14, with and without energy dissipation. When energy dissipation is not considered in the model, the joint dissipative moment of Eq. (34) is null, as opposed to the case in which energy dissipation is taken into account. Compared with the results in which the range of motion methodology was not applied (Fig. 10), the amplitude of motion of the center of mass of the main foot is significantly reduced in the anteroposterior direction ($x$-direction) for talocrural plantarflexion. This is true for both cases of the application of the range of motion methodology, that is, with and without energy dissipation, and for both ankle joint models (Fig. 14a-14d). For the superoinferior direction ($z$-direction), the same conclusion is reached (Fig. 14a, Fig. 14b, Fig. 14e, Fig. 14f). For the mediolateral direction ($y$-direction), concerning only the results without energy dissipation, the amplitude of motion for the proposed joint model is higher in the lateral direction (negative $y$) when compared with the case in which the range of motion methodology was not applied (Fig. 10b, Fig. 10c and Fig. 14c, Fig. 14e). This is also true for the Anderson and Pandy’s model [22], but the amplitude of motion is higher for both the medial and lateral directions (positive and negative $y$, respectively). The increase in the amplitude of motion occurs because, in the analysis in which the range of motion methodology was not considered, the center of mass of the main foot never exceeded the maximum latitude in the mediolateral direction ($y$-direction), therefore not violating the range of motion for inversion or eversion. When the range of motion methodology (described in Sect. 4) was introduced in the model, if an unacceptable position of the main foot was verified in any direction, a joint resistance moment was applied to place the main foot in an acceptable position. The direction and magnitude of this moment may have caused the main foot to reach an unacceptable position in the mediolateral direction, therefore reaching the maximum latitude for inversion or eversion. It is important to note that, for both ankle models and for the case of the range of motion methodology without energy dissipation, the direction of movement of the center of mass of the main foot in the anteroposterior direction ($x$-direction) (Fig. 14a) is similar to the results presented in Fig. 10a. However, for the superoinferior ($z$-direction) and mediolateral ($y$-direction) directions, the movement is non-regular due to the successive moment applications to prevent range of motion violation (Fig. 10b, Fig. 10c and Fig. 14c, Fig. 14e). In the results with energy dissipation, the amplitude of motion in the mediolateral direction ($y$-direction) is reduced when compared with the results without the range of motion methodology (Fig. 10b, Fig. 10c and Fig. 14d, Fig. 14f). With energy dissipation, it is clear the maximum latitude for plantarflexion is reached in the sagittal and transverse planes on both ankle models (Fig. 14b, Fig. 14d), and the center of mass of the main foot tends to move to the medial direction (positive $y$) (Fig. 14d, Fig. 14e). This is not consistent with the conclusions reached in Fig. 10 for the Anderson and Pandy’s model [22], where the center of mass of the main foot moved to the lateral direction (negative $y$). The direction and magnitude of the moment applied to restrict the plantarflexion range of motion may be responsible for this difference. In order to provide a better comparison between the two ankle joint models when the range of motion methodology was applied with and without energy dissipation, the snapshots for the first second of the simulation are presented in Fig. 15 and Fig. 16, respectively.

Fig. 14

Trajectory of the center of mass of the main foot performed by the two ankle joint models in the (a)-(b) sagittal, (c)-(d) transverse, and (e)-(f) frontal planes. Figures (a), (c) and (e) represent the results without energy dissipation, whilst (b), (d) and (f) consider energy dissipation. The green marker represents the initial position of the center of mass of the main foot segment. The $x$, $y$, and $z$ directions correspond to the anteroposterior, mediolateral, and superoinferior directions in the three cardinal planes, respectively (Color figure online)

Fig. 15

Snapshots of the movement of the biomechanical model for the first second of simulation using (a) the Anderson and Pandy’s [22] classical universal joint (See Online Resource 3 – 11044_2023_9955_MOESM3_ESM) and (b) the proposed modified universal joint (See Online Resource 4 – 11044_2023_9955_MOESM4_ESM). The range of motion methodology is considered without energy dissipation (Color figure online)

Fig. 16

Snapshots of the movement of the biomechanical model for the first second of simulation using (a) the Anderson and Pandy’s [22] classical universal joint (See Online Resource 5 – 11044_2023_9955_MOESM5_ESM) and (b) the proposed modified universal joint (See Online Resource 6 – 11044_2023_9955_MOESM6_ESM). The range of motion methodology is considered with energy dissipation (Color figure online)

The variation of the mechanical energy of both ankle models was also evaluated, and it is as depicted in Fig. 17. If the range of motion methodology is not considered, there is no energy dissipation, and, therefore, both systems are conservative. If the range of motion methodology is applied without energy dissipation, it can be noted that, at every joint resistance moment application, which is visible by the peaks present on the plots of the mechanical energy, the system loses energy. However, the system rapidly recovers since the energy is lost and restored in a short period of time. In this case, the system is conservative. If energy dissipation is considered, the systems dissipate energy. This analysis is valid for both ankle joint models. In addition, while the maximum latitude for plantarflexion is reached at around 0.73 seconds of simulation in the Anderson and Pandy’s model, in the proposed approach, the maximum latitude is reached earlier, at around 0.43 seconds. This phenomenon can be observed by the time at which the mechanical energy plateau is reached in the red plots of Fig. 17, and it means that the proposed approach dissipates energy faster than the Anderson and Pandy’s model.

Fig. 17

Variation of the mechanical energy for the (a) classical and (b) modified universal joints considered in the biomechanical model to model the ankle joint complex of the human foot (Color figure online)

The computational cost of the two ankle models analyzed in the biomechanical model was compared by assessing the number of function evaluations. Considering the cases in which the range of motion methodology was implemented with energy dissipation, when the 5-second simulations were used to evaluate the computational cost of the joint models, it was observed that the system showed numerical issues after it was fully dampened. The numerical issues are directly associated with the numerical error arising from the integration process that does not allow to obtain null velocities and, thus, fully stop the motion of the bodies. In turn, very low velocities are calculated, which provokes variations in the position and velocity of the bodies of the biomechanical model, leading to the application of a joint dissipative moment to the system. The dissipative moment causes substantial variations, in relative terms, of the coordinates of the biomechanical model, which reduce the time step required to meet the defined tolerances and proceed with the integration process, significantly decreasing the computational efficiency. Taking this information into consideration, Fig. 18 shows the computational cost for 1 second of simulation only, representing the simulation time before the occurrence of these numerical issues. Observing Fig. 18, it can be concluded that the proposed joint model is slightly less efficient when the range of motion methodology is not applied. Considering the situation in which the range of motion methodology is applied without energy dissipation, it can be concluded that the Anderson and Pandy’s model [22] is clearly more efficient than the proposed approach. The same tendency is observed when the range of motion methodology is applied with energy dissipation. In general, the introduction of energy dissipation to the system decreases the computational cost of both joint models. It is possible to conclude that the proposed approach does not present benefits in terms of computational cost, but it reproduces the anatomical location of the talocalcaneal articulation more realistically.

Fig. 18

Computational cost of the two joint models considered for the ankle joint complex of the human foot for 1 second of simulation time (Color figure online)

Finally, a third analysis was carried out, which aimed to determine the influence of the massless link length (variable $l$ on Eq. (20)) on the response of the biomechanical model. According to the study of Siegler et al. [80], which used computational tomography data of 26 healthy individuals aged between 18 and 35 years old to quantify the three-dimensional morphology of the talus, the height of this bone was found to be 0.0417±0.0057 m. Thus, in this study, four values for massless link length were tested, namely ($i$) 0.0000 m, which corresponds to the Anderson and Pandy’s model [22], ($ii$) 0.0360 m (0.0417−0.0057 m), (iii) 0.0417 m, representing the proposed approach, and ($iv$) 0.0474 m (0.0417+0.0057 m). The location of the talocrural joint was not altered, but the location of the talocalcaneal joint was adjusted to accommodate the variations in the massless link length. A schematic representation of these four cases in the sagittal plane is presented in the plots of Fig. 19.

Fig. 19

Schematic representation, in the sagittal plane, of the cases tested for the massless link length, namely (a) 0.0000, (b) 0.0360, (c) 0.0417, and (d) 0.0474 m. The talus bone is highlighted in grey (Color figure online)

Based on the first and second analyses (without and with the consideration of the range of motion methodology), it can be concluded that the use of the massless link to model the ankle joint complex of the human foot introduces substantial differences primarily on the mediolateral direction ($y$-direction) for all parameters studied. Thus, only the mediolateral direction was studied in the third analysis, and the obtained results for the center of mass of the main foot segment are depicted in the plots of Fig. 20. The range of motion methodology was not applied in this analysis. Observing Fig. 20, it can be concluded that, as the length of the massless link increases, the center of mass of the main foot segment moves progressively to the medial direction (positive $y$) on the transverse and frontal planes. The most distinct plot is the one in which the massless link length is equal to 0.000 m, corresponding to the Anderson and Pandy’s model [22]. In addition, although this phenomenon is not evidently identified by the observation of Fig. 20, the amplitude of motion of the center of mass of the main foot increases with the increase in the massless link length, as concluded by the successive increase in the width of the plots. With these conclusions in mind, it should be recognized that altering the initial conditions of the model or the parameters of the modified universal joint strongly influences the response of the model and leads to different results.

Fig. 20

Influence of the massless link length on the response of the biomechanical model. Trajectory of the center of mass of the main foot segment in the (a) transverse and (b) frontal planes. The $x$, $y$, and $z$ directions correspond to the anteroposterior, mediolateral, and superoinferior directions in the three cardinal planes, respectively (Color figure online)

7 Concluding remarks

A new skeletal model for the ankle joint complex of the human foot was presented in this work, the performance of which was examined and compared with the solution developed by Anderson and Pandy [22]. The main kinematic modeling aspects of these two joint approaches were described under the framework of multibody systems methodologies, taking into account the corresponding kinematic constraints, the resulting Jacobian matrix, and the right-hand side of the acceleration constraint equations vector. A methodology to restrict the range of motion in human articulations was explained and utilized. A demonstrative application example was considered using a biomechanical model of the right human foot and leg. The influence of the two joint modeling approaches on the response of the biomechanical model was studied. The two ankle models were compared using two different scenarios. In the first analysis, the biomechanical model was subjected to the gravitational force only, and in the second analysis, a methodology to restrict the range of motion in human articulations was utilized. An additional study was carried out to analyze the influence of the massless link length on the response of the biomechanical model. Overall, the introduction of the massless link resulted in substantial differences in the mediolateral direction for all variables studied, while the anteroposterior and superoinfoerior directions presented similar results. It can be concluded that, despite being less computationally efficient, the proposed joint model represents the anatomical location of the talocalcaneal articulation more realistically. Future studies should aim to validate the proposed modified universal joint using experimental data of the human movement. An appropriate marker set protocol based on the bony landmarks specified by the recommendations of the International Society of Biomechanics [87] must be developed, and the degrees-of-freedom of the considered biomechanical model must be guided utilizing the experimentally acquired data. The proposed joint model can be utilized to study both healthy and pathological populations. Since the main differences between the classical and modified universal joints are in the mediolateral direction, tasks that are strongly influenced by the inversion and eversion of the human foot, such as walking on inclined planes or performing three-point crutch-assisted gait, can be considered for validation.