Advecting tracer particles

    This code is just an example on how to have a velocity field moving along particles. A trouble with stability is that it is hard to realize when you get an eigenvalue, that this is a real flow with all its complexity and beauty. There are many modelisation steps to back-engineer. This is why I give a few examples of showing the actual velocity field of the eigenvectors, and it helps very much to let it advect particles, because we usually view the flow by the motion of particles instead of velocity vectors.

    clear all; clf; format compact
    % parameters
    n=50;       % gridpoints for the velocity field
    np=20000;   % numer of tracer particles
    Lx=2*pi;    % Domain size in x
    Ly=2*pi;    % Domain size in y
    dt=0.1;     % Time step
    tmax=30;   % final time
    % build the grid for the velocity field

    The particles

    We set many particles, located by their x and y coordinate. We chose this randomly. In addition to that, we will give a different color to the particles depending on wether they are on the right or on the left initially to help the eye see what happens and the global motion. For this selection, I use a logical array selp, which is true (value 1) if the particle is on the right and false (value 0) if it is on the left. We give the color blue to the “true” particles and the color red to the “false” ones. below in the plotting command.

    % Random initial positions of the particles
    % to select the particles
    % time loop
    for ind=1:length(tvec);

    The velocity field

    We build an artificial velocity field based on the streamfunction \displaystyle \phi=\cos(x)\cos(y) so that we have \displaystyle u=-\phi_y=\cos(x)\sin(y) and \displaystyle v=\phi_x=-\sin(x)\cos(y) and to let it be a little more interesting, we have the velocity field change in time instead as \displaystyle \phi=\cos(x-t)\cos(y) This velocity field is now defined on the cartesian grid that we have built before the loop using the very usefull function meshgrid; we now need to interpolate the field on this grid to know what is the velocity vector at the position of each particle. We do this using interp2 using the default linear interpolation scheme.

        % the velocity field
        % interpolate the velocity field


    To make things simple, we have a periodic box and a periodic velocity field and the particles may cross the boundaries of the box, so we need to find a way to insert them back on the other side of the box when they are advected out. This is made simply by using the function mod for the output of the very simple explicit advection scheme where the particle is advanced using the actual velocity for the time step dt


        % march one step
        % plotting
        axis equal; axis([0,Lx,0,Ly]);
        xlabel('x'); ylabel('y'); title('advection of particles');
        %print('-dpng','-r100',['frames/particles_' num2str(ind) '.png'])
    the last time of the simulation

    the last time of the simulation

    The code when it runs displays the figure synchronously so you see the animation, but please click on animation to see it saved for you:


    • Please put more than two colors
    • Please use another velocity field
    • Please build a velocity field using point vortices of a perfect fluid that advect each other in time and advect the particles using this instationnary field. You can as well have sources and sink, but then it will be nice to introduce new particles in the sources to avoid area without particles arround them.