Another important thing to note is that instabilities don't come from the constraints themselves, but from how constraint forces are distributed to the other constraints and how these constraints are geometrically related to each other. This is is the essence of the geometric stiffness concept: that constraint forces and constraint orientations can act as stiff forces (e.g. angular springs) in the transversal directions. Take as example a thread made of particles and massless rigid constraints between them: the constraint forces that act along constraints to restrict half of the degrees of freedom (DOFs), also act on the rest of the DOFs (i.e. the angles of the multi-pendulum) and couples them together (like in the chaotic behavior of the double pendulum).
Consider we are given a dynamical system with $n$ degrees of freedom and $m$ constraints. I will consider the case $m<n$ here, and leave the over-constrained case for future treatment. The constraint are given implicitly by $m$ constraint functions: $\psi(q)=0$. In fact, these functions measure the displacement along the constraint direction, so they can be considered as a coordinate.
The trick I will use is to consider a coordinate change from $q$ to another set of coordinates $(\psi, \theta)$. As you can see, we already know $m$ coordinates given by $\psi$, and we need to pull out the rest of $n-m$ coordinates $\theta$ from somewhere. This is quite simple to do for a thread (multi-pendulum) where we have closed formulas for all the coordinate changes, but not so straightforward for the general case. Turns out we don't need to know these nonlinear transformation formulas, but rather their Jacobians. This approach is inspired from an older way of solving constraints, known as the local coordinates approach. The main reference used here will be [Potra and Rheinboldt 1991] (On the Numerical Solution of Euler-Lagrange Equations).
The Jacobian of $\psi(q)$ is denoted by $J=\nabla_q \psi(q)$ and is usually easy to compute and necessary for constraint solving. This is not true though for $L$, the Jacobian of $\theta(q)$. But [Potra and Rheinboldt 1991] show that $L$ exists locally and can be computed numerically (through QR factorization of J). Moreover, it can be chosen so that $L^TL=1$ (i.e. orthonormal). Once we have $L$ (even if locally), we can start looking at the equations of motion and how they would look like in $\theta$-space (i.e. transversal direction). So let's see the equations of motion in maximal coordinates $q$ with only constraints and no other forces:
\begin{equation}
M\ddot{q} = J^T \lambda.
\end{equation}
Move the mass on the right-hand side and multiply to the left by $L$:
\begin{equation}
L\ddot{q} = L M^{-1} J^T \lambda.
\end{equation}
We do this because we know the fact that $L=\nabla_q \theta(q)$ and so $\dot{\theta}=L\dot{q}$. Differentiate another time and we get $\ddot{\theta}=L\ddot{q} + \dot{L}\dot{q}$. Therefore we get:
\begin{equation}
\ddot{\theta} = L M^{-1} J^T \lambda + \dot{L}\dot{q}.
\end{equation}
I will call the last term a gyroscopic term, as it only appears because of the nonlinear change of coordinates and I will ignore it from now on, because I know it's not the culprit for instability. But I do know that the term $L M^{-1} J^T \lambda$ is the one that's stiff. Most of this insight came (to me and others) from the first two papers mentioned in the first paragraph. So, in order to quantify the stiffness, we need to linearize this term (keep in mind we are already in a restricted vicinity of $q$ because of $L$). So let's differentiate the matrix $A = L M^{-1} J^T$ with respect to $\theta$:
\begin{equation}
\frac{\partial A}{\partial \theta}=\frac{\partial L}{\partial \theta} M^{-1} J^T + L M^{-1} \frac{\partial J^T}{\partial \theta}.
\end{equation}
But here's the catch: the transversal directions $L$ should not depend on the transversal coordinates $\theta$ but only on the other coordinates $\psi$. A similar choice of constraint directions $J$ can be made so that they only depend on $\theta$. What I cannot prove right now is that such a choice for both $J$ and $L$ always exists. However, I will assume it does, and continue based on that - namely, that $\frac{\partial L}{\partial \theta}=0$. Another thing we want to do, is express all terms as much as possible in terms of $\theta$, so we note that:
\begin{equation}
\frac{\partial J^T}{\partial \theta}= L^T \frac{\partial^2 \psi}{\partial \theta^2}.
\end{equation}
We can now finish our linearization, and after dumping the gyroscopic term, we get the following approximate local equation:
\begin{equation}
\ddot{\theta} = C + LM^{-1}L^T \left( \frac{\partial^2 \psi}{\partial \theta^2}\lambda \right) \theta = C + EK\theta,
\end{equation}
where $C$ is a constant term, $E$ is the inverse effective mass in transversal space and $K$ is the stiffness matrix as seen by the free degrees of freedom $\theta$. What we would like to show now is that there is a relation between this matrix $K=\frac{\partial^2 \psi}{\partial \theta^2}\lambda$ and the geometric stiffness $\tilde K = \frac{\partial^2 \psi}{\partial q^2}\lambda$. Given that $\frac{\partial^2 \psi}{\partial \psi^2} = \frac{\partial^2 \psi}{\partial \theta \partial \psi} = 0$, I was able to show that $\tilde K = L^T K L$. I may have to go again through the calculations, but it makes sense that the matrix $K$ is "pushed forward" in $q$-space and vice-versa.
Now, let's try and bring the last equation from transversal space back to $q$-space. Note that we linearized and approximated the original equations, so what we are bringing back now is basically just a portion of all the forces, but most importantly they are the stiff forces. This why we also leave out the constant term. Therefore, we multiply the equation to the left by $L^T$:
\begin{equation}
L^T\ddot{\theta} = L^T LM^{-1}L^T K \theta.
\end{equation}
We can now use the fact that $L^TL=1$:
\begin{equation}
\ddot{q} = M^{-1}L^T K \theta(q).
\end{equation}
We need again to linearize this equation (around a point $q^*$), so that stiffness shows its face:
\begin{equation}
\ddot{q} = M^{-1}L^T K L (q - q^*) + M^{-1}f^*.
\end{equation}
It is clear now to see that:
\begin{equation}
M\ddot{q} = \tilde K (q - q^*) + f^*.
\end{equation}
Again, we should ignore the constant terms and notice the fact that the stiffness matrix of the transverse forces in $q$-space is just $\tilde K$ - the geometric stiffness. Of course, this is just a linearization of the real force, but it is what matters locally, especially when doing numerical integration. Linear stability analysis is all we have, so the values of the $\tilde K$ matrix are very important for stability analysis in $q$-space (as the values of $K$ are in $\theta$-space). If the highest eigenvalue $\omega^2$ of $M^{-1}\tilde K$ does not satisfy $\omega < 2/h$, then the simulation becomes unstable. This is an idea that has already been proposed in [Andrews et al. 2017], but I was just not sure if it was true. Now that I've shown it is, it's time to focus on the modal analysis of $M^{-1}\tilde K$ as the proper way of studying stability in constrained based dynamics.
No comments:
Post a Comment