1.Define the simulation parameters

2.Initialize the simulation:

A.Particle locations: place NxN particles randomly on a two dimensional plane. E.g., determine the (x0, y0) coordinates for each particle using the rand() function. (the upper script 0 denotes the temporal index)

B.Particle velocities: zero for both velocity components. E.g., assign zero to (up0, vp0) for each particle suing the zeros() function.

3.Advance in time

for n=1, totaltimesteps

Compute the flow velocity at the corresponding location for each particle, ufn-1=cos(xn-1)sin(yn-1), vfn-1=-sin(xn-1)cos(yn-1)

Compute forces on each particle, fxn-1=6??f??R(ufn-1upn-1), fyn-1=6??f??R(vfn-1vpn-1)

Compute the particle velocity at n, upn=upn-1+?tfxn-1/(?pVp), vpn=vpn-1+?tfyn-1/(?pVp)

Compute the particle location at n, xpn=xpn-1+?tupn, ypn=ypn-1+?tvpn

end