Constraint force that keeps a body in a circular path with fixed radius

Question

I'm trying to simulate a simple electric motor from scratch. In order to do this I need to be able to apply forces to a body which can only travel around a circular parth, a rotating magnet. If my research is correct I now need something called constraints. I found this article which in the beginning describes pretty much exactly what I want to do:

https://www.toptal.com/game/video-game-physics-part-iii-constrained-rigid-body-simulation

However I am quite confused by the formular for the "Lagrange multiplier" as it doesnt enforce any kind of radius. I also tried implementing the constraint force like the article describes but the movement of my magnet now just seems random to me.

Please let me know what I'm missing. If you couldn't already tell by now I'm pretty new to this stuff and there very well could be a much easier solution to my problem.

Edit:

This is my magnet which is supposed to only be able to rotate around the origin. I've found two ways to make this possible.

1. To restrict the DOF of it and to only calculate the rotational acceleration produced by the acting forces. This is of course much easier and I think kind of hacky which is not really what I want.

2. I introduce some kind of constraint force which keeps the position of the magnet fixed.

During my simulation loop I'm adding all the acting forces to my acceleration (a = F/m) to then numerically integrate in order to estimate the next position of my object. I had hoped this constraint force would just be another force I'd be able to add to my acceleration.

Answer 1

$\def \b {\mathbf}$

General approach

Newton Equations of Motion \begin{align*} & m\,\ddot{\mathbf{R}}=\b F_a+\b F_c\tag 1 \end{align*} and the constraint equation e.g. for circular path \begin{align*} & x^2+y^2=r^2\quad,2\,x\,\dot{x}+2\,y\dot{y}=0\quad, \underbrace{\begin{bmatrix} x & y \\ \end{bmatrix}}_{\b C_c}\begin{bmatrix} \dot{x} \\ \dot{y} \\ \end{bmatrix}=0 \end{align*} with $~\b F_c=\b C_c^T\,\lambda~$ equation (1) \begin{align*} & m\,\ddot{\mathbf{R}}=\b F_a+\b C_c^T\,\lambda~\tag 2 \end{align*} choose $~x~$ to be the generalized coordinate , thus from the constraint eqaution you obtain

\begin{align*} &y=\sqrt{r^2-x^2}\quad\Rightarrow \dot{y}=-\frac{x}{y}\dot{x}\\\\ &\dot{\mathbf{R}}=\begin{bmatrix} \dot{x} \\ \dot{y} \\ \end{bmatrix}= \underbrace{\begin{bmatrix} 1 \\ -\frac{x}{y} \\ \end{bmatrix}}_{\b J}\,\dot{x}\quad \Rightarrow \ddot{\mathbf{R}}=\b J\,\ddot{x}+\frac{\partial \dot{\mathbf{R}}}{\partial x}\,\dot{x}= \b J\,\ddot{x}+\frac{\partial \left(\b J\,\dot{x}\right)}{\partial x}\,\dot{x} \end{align*} thus equation (2) \begin{align*} & m\,\b J\,\ddot{x}=\b F_a-\underbrace{m\,\frac{\partial \left(\b J\,\dot{x}\right)}{\partial x}\,\dot{x}}_{\b F_z}+\b C_c^T\,\lambda~\tag 3 \end{align*} multiply Eq. (3) from the left with $~\b C_c~$ you obtain \begin{align*} & \left[\b C_c\,\b C_c^T\right]\,\lambda=\b C_c\left(\b F_z-\b F_a\right)\tag 4 \end{align*} from here you get the generalized force $~\lambda~$ , and the constraint forces $~\b F_c=\b C_c^T\,\lambda~$

to obtain the equation of motions , multiply Eq. (3) from the left with $~\b J^T~$ \begin{align*} & m\,\b J^T\,\b J\,\ddot{x}=\b J^T\,\left(\b F_a-\b F_z\right)\tag 5 \end{align*}

notice that $~\b C_c\,\b J=\b 0$

• $~\b R~$ Position vector
• $~\b F_a~$ Applied force
• $~\b F_c~$ Constraint force
• $~\b \lambda~$ Generalized constraint force

---

Implementation

form the time derivative of constraint equation you obtain the constraint matrix $~\b C_c~$.

from the choose of the generalized coordinate $~(x~)~$ and the velocity equation $~\dot{\b{R}}~$ you obtain the Jacobian-Matrix $~\b J~$.

from here go to polar coordinate \begin{align*} &x=r\,\cos(\varphi)\quad,\dot{x}=-r\,\sin(\varphi)\\ &y=r\,\sin(\varphi)\quad,\dot{y}=+r\,\cos(\varphi)\\ \end{align*}

your generalized coordinate is now $~\varphi~$ thus $~\b C_c=\b C_c(\varphi)~$ and $\b J=\b J(\varphi)~$ with \begin{align*} &\b F_z=m\,\frac{\partial \left(\b J\,\dot{\varphi}\right)}{\partial \varphi}\,\dot{\varphi} \end{align*}

you can solve Eq. (4) for $~\lambda~$ and obtain the constraint force $\b F_c=\b C_c^T\,\lambda$

\begin{align*} \b F_c=\left[ \begin {array}{c} - \left( \cos \left( \varphi \right) {\it F_x}+\sin \left( \varphi \right) {\it Fy}+m\,r{\dot\varphi }^{2} \right) \cos \left( \varphi \right) \\ - \left( \cos \left( \varphi \right) {\it F_x}+\sin \left( \varphi \right) {\it F_y }+m\,r{\dot\varphi }^{2} \right) \sin \left( \varphi \right) \end {array} \right]\\ \end{align*}

with Eq. (5) you obtain the equations of motion \begin{align*} &m\,r^2\,\ddot\varphi=r\,(F_y\cos(\varphi)-F_y\sin(\varphi)) \end{align*}

The constraint force components $~(x,y)~$ system \begin{align*} &\b F_c= \left[ \begin {array}{c} -{\frac { \left( {\it F_x}\,x+{\it F_y}\,y+m{{ \dot x}}^{2}+m{{\dot y}}^{2} \right) x}{{r}^{2}}}\\ -{\frac { \left( {\it F_x}\,x+{\it F_y}\,y+m{{\dot x}}^{2}+m{{\dot y}}^{ 2} \right) y}{{r}^{2}}}\end {array} \right] \end{align*}

and the EQM's

\begin{align*} &\begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \end{bmatrix}=\b F_c \end{align*}

the initial conditions must fulfilled the constraint equation $~x_0^2+y_0^2=r^2~$

The simulation with $~r=1~,F_x=F_y=1,x_0=0.3~,y_0=\sqrt{1-0.3^2}$

Answer 2

There is an alternative way to the excellent answer @Eli provided.

Instead of using $(x,y,z)$ coordinates to track the position of each magnet, write the position in terms of coordinates that describe the motion. In this case of each magnet is constrained in the xy plane, but at a fixed radius $r$ from the center then you write

$$ \text{position vector} \\ \boldsymbol{p}(\theta) = \pmatrix{ r \cos \theta \\ r \sin \theta \\ 0} $$

where your position vector $\boldsymbol{p}$ is described in terms of ONE variable only, the angle $\theta$. This is called a generalized coordinate.

Now apply the derivative chain rule with fixed $r$ and varying $\theta$ to get velocity and acceleration

$$ \text{velocity vector} \\ \boldsymbol{v}(\theta, \dot{\theta}) = \pmatrix{ -r \dot{\theta} \sin \theta \\ r \dot{\theta} \cos \theta \\ 0} $$

$$ \text{acceleration vector} \\ \boldsymbol{a}(\theta, \dot{\theta}, \ddot{\theta}) = \pmatrix{ -r \ddot{\theta} \sin\theta - r \dot{\theta}^2 \cos\theta \\ r \ddot{\theta}\cos\theta - r \dot{\theta}^2 \sin \theta \\ 0 } $$

Now you decompose the velocity vector into two parts, based on your generalized coordinate speed, as so

$$ \boldsymbol{v} = \boldsymbol{J}\, \dot{\theta} = \pmatrix{ -r \sin \theta \\ r \cos \theta \\ 0} \dot{\theta} $$

where $\boldsymbol{J}$ is the jacobian which describes how much each velocity component changes with each change in the generalized coordinate.

Finally, a constraint force $\boldsymbol{F}_c$ must do no work in the system, so the power at any instant should be zero, $P = \boldsymbol{F}_c \cdot \boldsymbol{v} = 0$. This sets up the constraint equations for the forces

$$ \boldsymbol{F}_c \cdot \boldsymbol{J} = 0 $$

where $\cdot$ is the dot product, and in linear algebra world it is equal to $\boldsymbol{J}^\top \boldsymbol{F}_c =0$ where ${}^\top$ is the transpose operator.

So your system of equations are

$$\begin{gathered} \boldsymbol{F}_{\rm ext} + \boldsymbol{F}_c = m \boldsymbol{a} \\ \boldsymbol{J}^\top \boldsymbol{F}_c = 0 \end{gathered} $$

where $\boldsymbol{F}_{\rm ext}$ is any applied external forces.

In terms of components the above is

$$ \begin{aligned} Fx_c + Fx_{\rm ext} & = - m r \dot{\theta}^2 \cos \theta - m r \ddot{\theta} \sin \theta \\ Fy_c + Fy_{\rm ext} & = - m r \dot{\theta}^2 \sin \theta + m r \ddot{\theta} \cos \theta \\ Fz_c + Fz_{\rm ext} & = 0 \\ Fy_c r \cos \theta - Fx_c r \sin \theta & = 0 \end{aligned} $$

The system is to be solved in terms of the constraint forces $Fx_c$, $Fy_c$, $Fz_c$ and motion $\ddot{\theta}$.

The above process can easily be generalized for multiple bodies with multiple degrees of freedom, all formulated from your choice of generalized coordinates to track all the centers of masses.

$$ \boldsymbol{q} = \pmatrix{ \vdots \\ \theta_i \\ \vdots} $$

$$ \boldsymbol{v} = \pmatrix{ \vdots \\ \boldsymbol{v}_i \\ \vdots} = \mathbf{J}\, \boldsymbol{ \dot{q} } $$

$$ \boldsymbol{F}_c = \pmatrix{ \vdots \\ {\boldsymbol{F}_c}_i \\ \vdots} $$

$$ \boldsymbol{a} = \mathbf{J}\,\ddot{q} + \mathbf{\dot{J}} \dot{q} $$ $$ \mathbf{M} = \begin{bmatrix} \ddots & & \\ & m_i & \\ & & \ddots \end{bmatrix}$$ $$ \boldsymbol{F}_{\rm ext} + \boldsymbol{F}_c = \mathbf{M}\,\boldsymbol{a}$$ $$ \mathbf{J}^\top \boldsymbol{F}_c = \boldsymbol{0}$$