Application Center - Maplesoft

App Preview:

Flux of magnetic fields and gas expansion

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

surfintvec.mws

Surface Integrals of Vector Fields

Using Maple and the vec_calc Package

This worksheet shows how to compute surface integrals of vector fields using Maple and the vec_calc package. As examples we compute

* The Flux of a Magnetic Field

* The Expansion of a Gas

To start the vec_calc package, execute the following commands:

> restart;

> libname:="C:/mylib/vec_calc7", libname:

> with(vec_calc): vc_aliases:

> with(linalg):with(student):with(plots):

Warning, the protected names norm and trace have been redefined and unprotected

Warning, the name changecoords has been redefined

>

The Flux of a Magnetic Field

The magnetic field of a current i moving along an infinitely long straight wire has magnitude B = 1/4*mu[0]*i/(Pi*r) where mu[0] is a constant and r is the distance from the wire. Its direction is counterclockwise around the wire as given by the right hand rule. Suppose two wires are parallel to the z -axis and located at [x, y] = [2*a, 0] and [x, y] = [-2*a, 0] and each carries a current i but in opposite directions. We want to find the flux of the magnetic field through the rectangle -b <= x ` ` <= b , y = 0 , and 0 <= z ` ` <= c .

The magnetic field of the wire at [x, y] = [2*a, 0] is

> B1:=MF( [x,y,z], evall(mu[0]*i/(4*Pi*((x-2*a)^2+(y)^2))*[-y,x-2*a,0]) );

B1 := [proc (x, y, z) options operator, arrow; -1/4...

The magnetic field of the wire at [x, y] = [-2*a, 0] is

> B2:=MF( [x,y,z], evall(-mu[0]*i/(4*Pi*((x+2*a)^2+(y)^2))*[-y,x+2*a,0]) );

B2 := [proc (x, y, z) options operator, arrow; 1/4*...

The total field is

> B3:=MF( [x,y,z], evall(B1(x,y,z)+B2(x,y,z)) );

B3 := [proc (x, y, z) options operator, arrow; 1/4*...

For the purpose of plotting the functions (and making numerical computations) we make the following assignments to the variables:

> a:=1; b:=1; c:=1; i:=1; mu[0]:=1;

a := 1

b := 1

c := 1

i := 1

mu[0] := 1

We first plot the two functions separately:

> B1plot:=fieldplot3d(B1,-3..3,-1..1,-1..1, grid=[9,9,2], color=blue, arrows=SLIM, orientation=[-90,0], axes=normal, tickmarks=[7,3,3]):

> B2plot:=fieldplot3d(B2,-3..3,-1..1,-1..1, grid=[9,9,2], color=red, arrows=SLIM, orientation=[-90,0], axes=normal, tickmarks=[7,3,3]):

> display(array([B1plot,B2plot]));

[Maple Plot]

Notice that for both fields, the arrows point in the negative y -direction between the two wires. So the fields add up. Now we plot the total field:

> B3plot:=fieldplot3d(B3,-3..3,-1..1,-1..1, grid=[9,9,2], color=black, arrows=SLIM, orientation=[-90,0], axes=normal, tickmarks=[7,3,3]): B3plot;

[Maple Plot]

Before we go on, we unassign the parameters:

> a:='a'; b:='b'; c:='c'; i:='i'; mu[0]:='mu[0]';

a := 'a'

b := 'b'

c := 'c'

i := 'i'

mu[0] := 'mu[0]'

The rectangle may be parametrized by

> Rect:=[u,0,v];

Rect := [u, 0, v]

where -b <= u ` ` <= b , and 0 <= v ` ` <= c . The tangent vectors and the normal vector are

> Ru:=diff(Rect,u); Rv:=diff(Rect,v); N:=Ru &x Rv;

Ru := [1, 0, 0]

Rv := [0, 0, 1]

N := [0, -1, 0]

Notice that the normal points in the negative y -direction, which is the direction we saw in the plot that the magnetic field points.

On the rectangle, the magnetic field is

> B3(u,0,v); simplify(%);

[0, -1/4*mu[0]*i/Pi/(u+2*a)+1/4*mu[0]*i/Pi/(u-2*a),...

[0, mu[0]*i*a/Pi/(u^2-4*a^2), 0]

So the flux is

> interface(showassumed=0); assume(0<b,b<2*a,0<c);

> Muint(B3(u,0,v) &. N,u=-b..b,v=0..c); Flux:=value(%);

Int(Int(mu[0]*i*a/Pi/(-u^2+4*a^2),u = -b .. b),v = ...

Flux := -1/2*mu[0]*i*(ln(-b+2*a)-ln(b+2*a))/Pi*c

> Flux:=combine(%,ln);

Flux := mu[0]*i/Pi*c*ln(1/(sqrt((-b+2*a)/(b+2*a))))...

Note: There is no guarantee that if you re-execute these commands that Maple will give the answer in the same form.

Another way to compute this surface integral is to use the Surface_int_vector command (or its alias Siv ) from the vec_calc package which works directly with the parametrized surface and the vector field:

> Ruv:=MF([u,v],Rect);

Ruv := [proc (u, v) options operator, arrow; u end ...

> Siv(B3,Ruv,u=-b..b,v=0..c); Flux:=value(%);

Int(Int(mu[0]*i*a/Pi/(-u^2+4*a^2),u = -b .. b),v = ...

Flux := -1/2*mu[0]*i*(ln(-b+2*a)-ln(b+2*a))/Pi*c

There is a second way to compute this flux. The magnetic field is solenoidal because it is divergence-free:

> DIV(B3);

0

So the magnetic field has a vector potential, A .

> VEC_POT(B3,'A');

true

> A(x,y,z);

[-1/4*mu[0]*i/Pi/((x+2*a)^2+y^2)*(x+2*a)*z+1/4*mu[0...

By Stokes' Theorem, the flux is the line integral of the vector potential around the boundary curve. Here the boundary consists of four line segments:

> r1:=MF([u],[u,0,0]); #for u=-b..b

r1 := [proc (u) options operator, arrow; u end proc...

> r2:=MF([v],[b,0,v]); #for v=0..c

r2 := [proc (v) options operator, arrow; b end proc...

> r3:=MF([u],[-u,0,c]); #for u=-b..b

r3 := [proc (u) options operator, arrow; -u end pro...

> r4:=MF([v],[-b,0,c-v]); #for v=0..c

r4 := [proc (v) options operator, arrow; -b end pro...

So the flux is the total line integral:

> Liv(A,r1,u=-b..b)+Liv(A,r2,v=0..c)+Liv(A,r3,u=-b..b)+Liv(A,r4,v=0..c); Flux:=value(%);

Int(0,u = -b .. b)+2*Int(0,v = 0 .. c)+Int(mu[0]*i*...

Flux := 1/2*mu[0]*i*c*(-ln(-b+2*a)+ln(b+2*a))/Pi

>

The Expansion of a Gas

The velocity field of a gas is V = [x*y^2, y*z^2, z*x^2] . We want to find the expansion of the fluid out through the sphere x^2+y^2+z^2 = 4 which is defined as Expansion = Int(`V.`,S) ` ` = Int(Int(`V.N`,theta = 0 .. 2*Pi),phi = 0 .. P...

We enter the fluid velocity as a function:

> V:=MF([x,y,z],[x*y^2,y*z^2,z*x^2]);

V := [proc (x, y, z) options operator, arrow; x*y^2...

The sphere may be parametrized as

> Rsph:=[2*sin(phi)*cos(theta),2*sin(phi)*sin(theta),2*cos(phi)];

Rsph := [2*sin(phi)*cos(theta), 2*sin(phi)*sin(thet...

So the tangent vectors are:

> R[theta]:=diff(Rsph,theta);

R[theta] := [-2*sin(phi)*sin(theta), 2*sin(phi)*cos...

> R[phi]:=diff(Rsph,phi);

R[phi] := [2*cos(phi)*cos(theta), 2*cos(phi)*sin(th...

and the normal vector is

> N:= R[theta] &x R[phi];

N := [-4*sin(phi)^2*cos(theta), -4*sin(phi)^2*sin(t...

We simplify the third component:

> N[3]:=simplify(N[3]);

N[3] := -4*sin(phi)*cos(phi)

The components of N are negative in the first octant. So N points inward. We want the flux outward. So we reverse N :

> N:=-N;

N := [4*sin(phi)^2*cos(theta), 4*sin(phi)^2*sin(the...

On the sphere, the velocity is

> VR:=V(op(Rsph));

VR := [8*sin(phi)^3*cos(theta)*sin(theta)^2, 8*sin(...

Note: We need the op command to strip the square brackets off Rsph .

The dot product of the velocity field and the normal is:

> VN:=VR &. N;

VN := -32*(-1+cos(phi)^2)*sin(phi)*(cos(theta)^2-co...

So the expansion is:

> Muint(VN,theta=0..2*Pi,phi=0..Pi); Expansion:=value(%);

Int(Int(-32*(-1+cos(phi)^2)*sin(phi)*(cos(theta)^2-...

Expansion := 128/5*Pi

Another way to compute this surface integral is to use the Surface_int_vector command (or its alias Siv ) from the vec_calc package which works directly with the parametrized surface and the vector field:

> R:=MF([theta,phi],Rsph);

R := [proc (theta, phi) options operator, arrow; 2*...

> Siv(V,R,theta=0..2*Pi,phi=0..Pi); Expansion:=value(%);

Int(Int(32*sin(phi)*(-cos(theta)^2+cos(theta)^4+2*c...
Int(Int(32*sin(phi)*(-cos(theta)^2+cos(theta)^4+2*c...

Expansion := -128/5*Pi

Notice there is a minus sign difference in the answer because we did not reverse the normal.

There is a second way to compute this expansion. By Gauss' theorem, the expansion is the volume integral of the divergence of the velocity over the interior of the sphere. In this case, the divergence is:

> DIV(V);

proc (x, y, z) options operator, arrow; y^2+z^2+x^2...

which in spherical coordinates is

> divV:=rho^2;

divV := rho^2

The spherical volume Jacobian is

> J:=rho^2*sin(phi);

J := rho^2*sin(phi)

So the expansion is

> Muint(divV*J,rho=0..2,theta=0..2*Pi,phi=0..Pi); Expansion:=value(%);

Int(Int(Int(rho^4*sin(phi),rho = 0 .. 2),theta = 0 ...

Expansion := 128/5*Pi

>

>