Streamfunction of a vector field

General Discussion
1

Hello,

I just discovered Freefem and use it for my master degree in bioscience.
I’d like to visualize the streamlines of a 2d vector field so that their isolines must be always orthogonal to the vector field. The “streamfunction” S satisfied the “vorticity” equation : div(grad(S))=-curl(F)[z] since grad(S)=[-Fy,Fx].

Unfortunalty i can’t get streamlines orthogonal to the isolines, but not so far, do I miss something ?
My example is very simple, it’s a square with three holes. First I solve the laplacian equation using :

  • dirichlet bc on borders 1,2,3,4,5
  • zero flux bc on border 6
    Then I solve the streamfunction :
  • vort is the z component of the curl(F)
  • borderDerivative is the normal derivative using the definition grad(S) = [-Fy,Fx]

Thanks for the time to help.

Hélène Ségalier.

PS : I’m not a mathematician (a biologist who learn math sometimes)

Here my edp file

load "gmsh"
load "iovtk"

mesh Th = gmshload("square.msh");

fespace Ph(Th, P2);

int[int] reg = regions(Th);
int[int] tags = labels(Th);

cout<<reg<<endl;
cout<<tags<<endl;

Ph u,v;

problem laplacian(u, v,solver=UMFPACK,eps=1e-5)
= int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))
+ on(1, u=0)
+ on(2, u=0)
+ on(3, u=1)
+ on(4, u=1)
+ on(5, u=1)
;
laplacian;

Ph Fx,Fy;
Fx = -dx(u);
Fy = -dy(u);

Ph S;
func rotZ = -dx(Fy) + dy(Fx);
func borderDerivative = (-Fy*N.x + Fx*N.y);

problem stream(S,v,solver=UMFPACK,eps=1e-5)
= int2d(Th)(dx(S) * dx(v) + dy(S) * dy(v))
- int2d(Th)(rotZ*v)
- int1d(Th)(borderDerivative*v);

stream;

real SMin = S[].min;
cout<<"S min = "<<SMin;
S = S-SMin;

real SMax = S[].max;
cout<<"S max = "<<SMax;
S = S/SMax;


int[int] Order = [1, 1];
string DataName = "u S";
savevtk("export.vtu", Th, u, S, dataname = DataName, order = Order, bin = false);
2

You can see in example
cavityNewton.edp (2.8 KB)

TH boundary condition on the `

// Problem stream-lines (with solve)
solve streamlines (psi, phi)
	= int2d(Th)(
		  dx(psi)*dx(phi)
		+ dy(psi)*dy(phi)
	)
	+ int2d(Th)(
		- phi*(dy(u1) - dx(u2))
	)
	+ on(1, 2, 3, 4, psi=psiBC);

where psiBC is the psiBC(x(s),y(s)) = \int_0^s U.n ds
where s is a curve abcisse on border and U is the flow field.

3

Hello @frederichecht,

I really want to thank you for the time you took to answer.

  • I’ve already showed this example cavityNewton, in that example the velocity field is normal to the borders (slip condition), so that the boundary integral is 0 and psi=0 on the boundaries.

  • In my case (a square with 3 holes) that integral is not obvious and i don’t know how to compute it with freefem. Moreover one must choose an arbitrary origin point on the border ?

  • I’m not sure to understand why you use a dirichlet bc instead of using the border integral in the weak formulation as follow :
    image
    Could you explain ?

Best regards.

4

1)Because in this case F.T = 0 so the integral is zero . (T == tangent )
2) the Dirichlet BC. are more precise , and you can see recirculation zone close to the corner.

So I make a mistake on the previous mail the psiBC(s) = int_0^1 F.T not int_0^1 F.N

Best Regards,

Frédéric Hecht.

5

-Yes F.T=0 !

  • What do you mean by recirculation ?
    I cant see how to compute psiBC for each node on the border to prescribe the dirichlet bc on(gamma, psi=psiBC)

Thanks for your support.

6

Of coarse, if you do not know the way to compute de Dirichlel BC you can use Neuman BC.