I am building a biomechanical model that consists of spherical joints. I wanted to represent these spherical joints by quaternions. I am struggling to derive how the muscles contribute to the dynamical system.
When you model a system where your joints are represented by Euler or Tait-Bryan angles, you take a derivative of the muscle length w.r.t. generalized coordinate and that is your moment arm ( $R_{q_i} = \partial l(\textbf{q})/\partial q_i$ ). It is based on the virtual work equivalence ($sdl = \tau_{\phi}d\phi$, where $s$ is muscle tension force, $l=l(\phi)$ is path length along which tension $s$ is acting). You just multiply it by the muscle force and you can add it to your EoMs as generalized forces like this:
$\textbf{M} \ddot q = \textit{f} + R*muscleForce $.
The problem is, how to apply this method (or another approach that would work) when using quaternion as an representation of rotation. The muscle length is now a function of quaternion and the equations of motion now have angular velocities instead of generalized speeds, so now the muscle should probably act as an external torque in the body's frame.
As quaternion basically represents the axis and angle of a rotation, so there may be a way how to calculate the axis around which the torque is acting based on the current length of the muscle, but I haven't found a way how to do it. Thanks for help.
Edit:
So I have found out the way how to do it. The tendon excursion method is basically calculating the jacobian. So I tried to calculate it the same way with the quaternion representation. The jacobian of the tendon length $J$ is now $J = \partial l(\textbf{Q})/ \partial Q_i$, where $\textbf{Q}$ is a quaternion and $Q_i$ are the elements of quaternion. Since I have my equations of motion in terms of angular velocities, I need to map the jacobian from generalized coordinates to spatial coordinates of the rigid body. I used the mapping based on this article - https://arxiv.org/abs/0811.2889:
$\omega' = 2G\dot{\textbf{Q}}$,
$G = \begin{bmatrix}-Q_1 & Q_0 & Q_3 & -Q_2\\ -Q_2&-Q_3& Q_0& Q_1 \\ -Q_3& Q_2& -Q_1& Q_0\end{bmatrix}$
where $\omega'$ is angular velocity in body reference frame, $\dot{\textbf{Q}}$ is time derivative of quaternion and $G$ is the mapping matrix. The correct spatial moments acting on the rigid body from the muscles forces are now:
$M^B = \frac{1}{2} G \frac{\partial l(\textbf{Q})}{\partial Q_i} * MuscleForce$
I have checked that the results with this method are correct. The problem is I don't actually know why this is now divided by 2 instead of multiplying by 2 as I have came to this fraction experimentally.
