Archive for the physics Category

About Coulombic Friction models

Posted in code, physics with tags , , on August 16, 2010 by maxpower3141

So I have been lately somewhat involved in modeling Coulombic friction in my (perhaps game-) physics simulations project, as you can see from my post on writing test code and the first O’Caml rant.

Just today I spent some time in augmenting the test-case for the static friction case, and it occurred to me, that you, my loyal reader and supporter, might be interested in what is involved with the math and code when implementing these kind of models, so here we go again:

First of all, let us set the background: We are interested in this case friction between two bodies connected with some kind of constraint – let this constraint be realized by the equation: \phi(x_i) = 0. The function \phi(x_i) could be almost anything involving the coordinates x_i, i = 1, 2, 3, \ldots – for example in my test-case has \phi(x, y) = \sin x + y = 0, which just tells that whatever the rest of the dynamics are, the point-body described by the coordinates (x, y) should always lie on the negative sine-curve.

Then basic Lagrangian mechanics tells us that there exists a quantity \lambda, which ensures that when one applies the force(s) F_i = \lambda \frac{\partial \phi}{\partial x_i} to the body (bodies), then the constraint can be satisfied. This force is the contact or normal force of the constraint and it has the special property that it doesn’t perform any physical work, because it is always orthogonal to the movement-vector of the body.

Then the typical, and quite accurate model, for friction force caused by the constraint \phi is that the movement of the body (/ the relative movement of the bodies) is slowed down by a force that is linearly proportional in magnitude to the normal force of the constraint F_i = \lambda \frac{\partial \phi}{\partial x_i}. Then of course the magnitude of the constraint force is \| F_i \| = | \lambda | \sqrt{\sum_i \left(\frac{\partial \phi}{\partial x_i}\right)^2} and the proportionality constant of the friction force is slightly higher in case the body is not (bodies are not) moving.

Until now, I have mostly concentrating on the dynamical friction, where the friction force is completely opposite to the velocity vector of the bodies v_i:
F^{fric}_i = - \xi \frac{v_i}{\|v\|} \| F_i \| = - \xi \frac{v_i}{\|v\|} | \lambda | \sqrt{\sum_i \left(\frac{\partial \phi}{\partial x_i}\right)^2}, where \xi is the coefficient of the friction.

Now we get to the actual predicament: How to do this in practice? The force is proportional to another force, but alas, it is the absolute value that it is proportional to! We would get perfect model for linear dependence – meaning something like F^{fric}_i = \lambda \ldots, but we do have linear dependency on the velocity – almost: it is linearly dependent on the direction of the velocity.

In order to study this, I devised three models how to numerically solve them:

1) Make the linear dependency on the velocity linear and use other values from previous known state – this results effectively in:
F^{fric}_i = - \xi \frac{v_i}{\|v_0\|} | \lambda_0 | \sqrt{\sum_i \left(\frac{\partial \phi}{\partial x_i^0}\right)^2}, where v_0 is the previous velocity, \lambda_0 is the previous Lagrange multiplier and x_i^0 are the coordinates from previous state. The nice thing about this model is that it is really simple and computationally inexpensive to implement, but it has the weird property that the velocity used comes from two different points in time, but according to my tests the model is quite ok (except for really smally velocities).

2) Make the force linearly dependent on the Lagrange multiplier \lambda, \rightarrow
F^{fric}_i = - \xi \frac{v_i^0}{\|v_0\|} \lambda \  sign (\lambda_0)\  \sqrt{\sum_i \left(\frac{\partial \phi}{\partial x_i^0}\right)^2}
Here of course |\lambda| = \lambda \  sign (\lambda_0) is always true, except when the constraint force changes directions between two time-steps, but in these situations normally anyways the bodies can drift apart (in case we have a normal contact-constraint \phi(x_i) \ge 0) and the model works nicely as it is completely up-to-date with the constraint force, but alas it is somewhat more complex to compute.

The third option, that I haven’t tried out yet, is to just throw away all quantities of current step and just throw the entire force on the right-hand-side of the dynamics equation, making it basically a constant force over each time-step and this just might be completely sufficient for my purposes but all this still requires some work and testing, as the dynamics model I’m using is somewhat novel. Hopefully I can share more about that later on. For now, I bid you good night.

Btw, I’m planning on modeling the rest-friction case as just conditional constraint in the tangent-space of the original constraint. More on that to come. 🙂


Writing test-code

Posted in code, physics, visualization with tags , on August 11, 2010 by maxpower3141

One of the tasks, that in my knowledge is too often neglected (well you know I do this!), is writing test-code for new functionality.

Recently in my post on an O’Caml rant, I was complaining how O’Caml made it impossible for me to do the “simplest of all computer operations” – sleep. I simply wanted to make my CPU(s) idle for a while after drawing some visualization graphics, in order to produce a fluid animation, and it seemed hopeless. I just went on and implemented (the dreaded) busy-loop. 🙂

This test-case, even though it provided visual confirmation that my dynamic friction model (well the first option anyways – I have 2 others to check still) was ok, but it was lacking something important: Error analysis – the Final Frontier for programmers, where even the uber-nerds wont go without some serious recommending.

There are a few easy things we can check out in these sort of cases:

1) Distance from constraint manifolds (for both position and velocity variables)
2) Constraint force (I actually projected the total force vector to constraint-manifold normal. subtracted gravity before projection in order to measure directly the force felt by particle due to constraint)
3) Friction force (should depend linearly on the constraint force)

Computing these values turned out to be already a task, but visualizing them over time turned out even larger task – here’s the gif-animation:

Animating particle sliding on the sine curve

Of course the classical “O’Caml list-reverse” happened, so that the plot is like a polygraph, lie-detector type of curve (For me this happens all the time with O’Caml lists – and of course any iteration through the list inverts it! :))

But furthermore the results seem good, except one tiny thing: The friction force seems to diverge when going through zero-velocity of the particle – how very odd, no? I spent already some time verifying the math – it all matches, so what’s wrong? I know my friction model in this case is not perfect, but the discrepancy should not be so noticeable.

Well, the (yet again too obvious) solution was that this model had only implemented dynamical friction!

    Rest friction not yet implemented!!

. D’Oh! 🙂 The classical – without rest friction the plot is supposed to look something like that!!!

So well, moving on to implement rest friction and then to try the two other friction models! I’ll be back with more cool things! … Like more code… And physics.. and math.. and.. Uh.. Yeah – I know – super cool!Well, I am having a glass of wine while coding now, so, I guess I’m not completely hopeless, yet? Ok? And if I play a tad on my guitar?


Nice particle-based physics engine

Posted in code, physics on August 4, 2010 by maxpower3141

Geeks3d has an article on Lagoa physics engine, which seems quite impressive.

Direct quote from site:

Lagoa Multiphysics is a new physics engine created by Thiago Costa. This physics engine is focused on particles.

I have been thinking about this kind of idea to approximate continuous matter as a set of constrained point-bodies (particles) and here mr. Costa brilliantly shows how impressive things you can create with this kind of system. Hopefully later on he will add more rigid bodies as well and destruction thingies.

Also this kind of solver should be just great for GPU-acceleration for physics!

Ocaml strikes back

Posted in code, physics, visualization on August 4, 2010 by maxpower3141

Grr – I got angry with my trusty old O’Caml’s back again.

It all started with my hobby physics engine project (I’ll talk about this more at some later post*) that I’ve been putting together in my free time – I wanted to visualize how well my first draft of joint-friction works; Should be easy with a powerful language, no sweat, right? Well almost. 🙂

*) Actually now that I think about it, I’ll probably never dedicate a whole post to explaining it, but you, my loyal reader, can gather bits and pieces of it from various future posts, or something. 🙂

I already had most of the dummy code for my test-case written (point-particle on a sine-curve)


type state = { x : float; y : float; vx : float; vy : float; lambda : float };;
let h = 0.02;;
let g = 10.0;;
let mass = 1.0;;
let fricCoef1 = 0.1;;

Then some update code (manully solving the constraint equation for phi(x,y) = y + sin x = 0 etc. and then the fun part:

Visualizing the slide!

Grab google and O’Caml reference manual ( and we have all we need:

(In case you are on Linux (more on that later 🙂 type something like ocaml graphics.cma to get the interpreter with X11 support)

Open Graphics module, create default window surface and resize it to (1024×768):

open Graphics;;
open_graph "";;
resize_window 1024 768;;

Then define starting point on window (lower right being 0,0!), with which factor to scale the gfx and draw the basic coordinate axles (x = 0 and y = 0)

let scale = 100.0;;
let x0 = 100;;
let y0 = 400;;

let drawAxis () =
set_color (rgb 0 0 0);
set_line_width 2;
moveto x0 10;
let w = size_x () in
let h = size_y () in
lineto x0 (h-11);
moveto 10 y0;
lineto (w -11) y0;

Easy, no?

Then just define the functions needed to draw the sine-curve, then the whole background and ultimate the whole state:

let pixelcoordX x = x0 + int_of_float (scale *. x);;
let pixelcoordY y = y0 + int_of_float (scale *. y);;

let drawCurve () =
let nPoints = 50 in
set_color (rgb 120 120 120);
set_line_width 1;
let first = ref true in
for i = 0 to nPoints do
let x = -.0.7 +. (float_of_int i) *.
9.6 /. (float_of_int nPoints) in
let y = -.(sin x) in
let drawCall =
if !first then
drawCall (pixelcoordX x) (pixelcoordY y);
first := false

(* drawCurve ();; *)

let drawBackGround () =

let drawState state =
let x = state.x in
let y = state.y in
set_color (rgb 220 20 120);
set_line_width 3;
draw_circle (pixelcoordX x) (pixelcoordY y) 6

Almost there! Only thing missing is some loop to render the frames – what this means, is that I need to just draw my graphics to the windowing-system backbuffer, post that to the frontbuffer and wait a while to get smooth animation – preferrably so that the the time used to compute physics step and draw the graphics is comparable to physics-stepsize.

But wait, there’s no sane way to wait in O’Caml!!! The Unix-module has wait-function for complete seconds only! So – I can’t wait 0.02 seconds which is my physics-stepsize. Unix-select could apparently also be used, but for some reason [] [] [] 0.02 causes an exception (The system-call is pre-empted by an interrupt apparently). So what is the solution? Apparently there isn’t any? How is one supposed to code this thing? I mean, these are preeeettyyy basic OS-services that should be always available. Well – after googling and trying out desperate things for an hour or so, I deviced an elegant solution: The busy-loop – here we go:

let drawFrame state =
auto_synchronize false;
(* ignore ( [] [] [] 0.01); *)
(* ignore (Unix.sleep 1); *)
for i = 1 to 90009 do (); done;
drawState state;
auto_synchronize true;

let runNFrames n state =
let s = ref state in
for i = 0 to n do
drawFrame !s;
s := fricstep !s

runNFrames 1000 s;;

Beautiful, isn’t it? Well, with this trickery I was able to play my smoooth animation and verify that visually the quality is ok and that the test-particle stops quite nicely to the bottom of the hill.