3
$\begingroup$

I am currently working on an impacting system which is basically just a spring damper and a circular enclosure. Because of the rotational symmetry of the problem I need the vector field in polar coordinates to derive the sticking flow along the enclosure. It should be an easy problem, but contrary to my expectation there are no scores of textbooks which contain the solution and I am not used to mechanical problems.

The system in Cartesian coordinates is the following: $$ \left( \begin{array}{c} \dot{x} \\ \ddot{x} \\ \dot{y} \\ \ddot{y} \\ \end{array} \right)= \left( \begin{array}{cccc} 0 & 1 & 0 & 0 \\ -\nu^2 & -\gamma & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -\nu^2 & -\gamma \\ \end{array} \right) \left( \begin{array}{c} x \\ \dot{x} \\ y \\ \dot{y} \\ \end{array} \right) $$

I have replaced $x$, $y$ and their derivatives by $x=r \cos(\theta)$, $y=r \sin(\theta)$ and their derivatives as was also done in this question. Solving for $\ddot{r}$ and $\ddot{\theta}$ I obtained $$ \ddot{r}= -r \nu ^2-\gamma \dot{r}+r \dot{\theta} ^2 $$ and $$ \ddot{\theta}=-\frac{\left(r \gamma +2 \dot{r}\right) \dot{\theta}}{r} $$ which I verified with mathematica.

But when I test using the following Matlab-code

nu2=1.96; gam=0.28;
ODE1=@(t,z) [0 1 0 0;-nu2 -gam 0 0;0 0 0 1;0 0 -nu2 -gam]*z;
ODE2=@(t,z) [z(2); 
            -nu2*z(1)-gam*z(2)+z(1)*z(4)^2; 
             z(4); 
            -z(4)*(gam*z(1)+2*z(2))/z(1)];
z0=[1;0;1;0.5];
tspan=[0, 1];
z0p=[1.4142 0.3536 0.7854 0.3536];
[T1,Z1]=ode45(ODE1,tspan,z0);
[T2,Z2]=ode45(ODE2,tspan,z0p);
Z2c=[Z2(:,1).*cos(Z2(:,3)),Z2(:,1).*sin(Z2(:,3))];
hold on
plot(Z1(:,1),Z1(:,3));
plot(Z2c(:,1),Z2c(:,2));
hold off
legend('Cartesian','Polar')

I find that they are not equivalent. Now I am wondering whether there is a conceptual error as I am fairly confident in my derivation. So I am starting to think that my approach must be wrong and I would happy if someone could tell me where my error lies.

Update

Found the problem, answer below.

$\endgroup$
2
  • $\begingroup$ What exactly do you mean by "they are not equivalent" ? The plots are different or what? Note that while the first system is very good for simulations (its right hand side is $C^{\infty}$ vector functions), the second one has singularity at $r = 0$ which (in my opinion) can spoil calculations a lot. $\endgroup$
    – Evgeny
    Jul 23, 2016 at 16:54
  • $\begingroup$ Yes, I meant different plots, but I found my error and updated accordingly. Thank you, I am aware of the singularity and initially suspected my problems somewhere with it. But in my setting $\dot{r}$ is always normal and $\dot{\theta}$ is always tangential to my impact boundary, which makes things very convenient. $\endgroup$
    – crown42
    Jul 23, 2016 at 18:32

1 Answer 1

Reset to default
0
$\begingroup$

It turns out the problem was in the calculation of the initial velocities which should be $$ \dot{r}=\dot{x}\cos(\theta)+\dot{y}\cos(\theta) \quad \text{and} \quad\dot{\theta}=\frac{\dot{y}\cos(\theta)-\dot{x}\cos(\theta)}{r} $$ Then the the system above is correct and both formulations yields the same plot of the correct initial conditions in polar coordinates are given which are

z0p=[1.4142    0.3536    0.7854    0.2500]
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .