Incompressible and irrotational velocity fields

In this worksheet, we

(1) compute the velocity potential of an irrotational fluid and graph its contour plot,

(2) compute the stream function of an incompressible fluid and graph its contour plot,

(3) plot the field plot of the velocity field of a fluid then solve the corresponding system of ODEs to determine the motion of a fluid particle,

(4) compute the circulation of an irrotational fluid about a closed curve (which is zero).

sm311o_13.mws,wdj,3-27-98

> restart;
with(linalg):
with(plots):

The velocity potential of an irrotational fluid

> v:=(x,y)->vector([x, -y]):
v(x,y);

This tells represents the velocity field of a fluid which (as we will see) is both irrotational and incompressible.

> v_new:=(x,y,z)->vector([x, -y,0]):
v_new(x,y,z);
curl(v_new(x,y,z),[x,y,z]);

This tells us that the velocity field is irrotational (or conservative). Another way to cheack that v(x,y) is irrotational is to check that dM/dy = dN/dx.

> fieldplot(v(x,y),x=-1..1,y=-1..1);

> potential([v(x,y)[1],v(x,y)[2]],[x,y],`phi`);
Is this velocity field irrotational?

> vel_potential:=(x,y)->phi:
vel_potential(x,y);

> grad(vel_potential(x,y),[x,y]);
Is the gradient of the velocity potential equal to the original velocity field?

>

> contourplot(vel_potential(x,y),x=-1..1,y=-1..1);

The stream function of an incompressible fluid

> v:=(x,y)->vector([x, -y]):
v(x,y);
diverge(v(x,y),[x,y]);

This tells us that the velocity field is incompressible.

> potential([-v(x,y)[2],v(x,y)[1]],[x,y],`psi`);
Is this velocity field incompressible? If true, let psi denote the stream function.

> stream:=(x,y)->psi:
stream(x,y);

> gradplot(stream(x,y),x=-1..1,y=-1..1);

> contourplot(stream(x,y),x=-1..1,y=-1..1);

The trajectory of a particle in the fluid

Now, let us simply take a particle at a point (x(t),y(t)) moving with velocity vector v(x,y)=<x,-y> and compute its trajectory:

> with(DEtools):

> sys_velocity:=diff(x(t),t)=x(t),diff(y(t),t)=-y(t);

> DEplot({sys_velocity},{x(t),y(t)},t=-2..2,x=-1..1,y=-1..1, arrows=SMALL,title=`system of ODEs for (x',y')=v(x,y)`,color=green);

These particles are moving along the contour lines of the stream function.

Now, let us simply take a particle at a point (x(t),y(t)) moving with velocity vector stream(x,y)=<y,x> (note this is orthogonal to f(x,y)) and compute its trajectory:

> sys_stream:=diff(x(t),t)=y(t),diff(y(t),t)=x(t);

> DEplot({sys_stream},{x(t),y(t)},t=-2..2,x=-1..1,y=-1..1,arrows=SMALL,title=`system of ODEs for (x',y')=grad(stream(x,y))`,color=blue);

These particles are moving along the contour lines of the velocity potential function.

The circulation of an irrotational fluid

This section requires a package to compute line integrals written by Joe Riel,
<jsr@sparc.sandiegoca.ncr.com>.

Finally, we shall see that the circulation of v(x,y) is 0 about a particular path C. This is consistent with the fact that v(x,y) is irrotational.

> read `c:/oldd/maplev4/share/analysis/pathint.mtx`;
This command reads in the file pathint.mtx (which you must have saved in the directory above or you must edit this path above to show the directory where you have loaded pathint.mtx). The file pathint.mtx will open on your screen. You don't need it so you may either minimize it or close it and then continue with the commands below.

>

> Gamma1:=[[t,2*t^2],t=0..1];
r:=[x,y];
v:=(x,y)->[x,-y];
Pathint(dot(v(x,y),D(r)),r=Gamma1);
pathint(dot(v(x,y),D(r)),r=Gamma1);

> Gamma2:=[[t,2*t],t=0..1];
Pathint(dot(v(x,y),D(r)),r=Gamma2);
pathint(dot(v(x,y),D(r)),r=Gamma2);

Let C be the closed curve represented by first traveling over Gamma1 then traveling over Gamma2 in the opposite direction (this is called -Gamma2). We write this as C = Gamma1 - Gamma2. Next, we compute the circulation of v about C.

The circulation of the velocity field v=< M,N > (which is equal to the flux of < -N,M >, by Green's theorem) about the curve C measures the net amount of the field < -N,M > passing through the curve C, regarded as a "screen". As the picture below shows, this appears to be zero since each arrow in the fieldplot for < -N,M > = < y,x > which is "captured" by the upper part of the screen (represented by Gamma) passed through the lower part (represented by Gamma), so the net amount captured is zero.

> P0:=fieldplot([y,x],x=-.5..1.5,y=-.5..2.5):
P1:=plot([t,2*t^2,t=0..1],color=red,thickness=2):
P2:=plot([t,2*t,t=0..1],color=green,thickness=2):
display([P0,P1,P2]);

Another example, and a picture of all the curves at once:

> Gamma3a:=[[t,0],t=0..1];
Gamma3b:=[[1,2*t],t=0..1];
Pathint(dot(v(x,y),D(r)),r=Gamma3a)+Pathint(dot(v(x,y),D(r)),r=Gamma3b);
pathint(dot(v(x,y),D(r)),r=Gamma3a)+pathint(dot(v(x,y),D(r)),r=Gamma3b);

> P1:=plot([t,2*t^2,t=0..1],color=red):
P2:=plot([t,2*t,t=0..1],color=green):
P3a:=plot([t,0,t=0..1],color=blue,thickness=3):
P3b:=plot([1,2*t,t=0..1],color=blue,thickness=3):
T1:=textplot([.3,1,`Gamma1`]):
T2:=textplot([.6,.35,`Gamma2`]):
T3:=textplot([1.05,.7,`Gamma3a+Gamma3b`],align={RIGHT}):
display([P1,T1,P2,T2,P3a,P3b,T3]);

>

>