Can pdetool() output the stream function?

Hello,

I am struggling to work out a matlab program for the numerical calculation of the streamfunction directly from the potentialfunction. I have the latter as a solution of a 2D Laplace equation from the pdetool() and would like to generate contour plots of both the potentialfunction as well as the streamfunction to form finally a net of coordinates that are rectangular in every point.

So far I tested flowfun.m from Kirill K. Pankratov, but this has problems when the function is not continuous. The 2D region where my PDE is defined has however several holes and an integration of the Cauchy-Riemann equations is numerically difficult at these holes.

Streamlines (that are easily provided by matlab) are not helpful because I have to specify starting points. What I need however is the scalar value of the streamfunction at every point of the 2D region.

Possibly the streamfunction (as the harmonic dual of the potentialfunction) can be derived directly during the FEM calculation of the solution of the PDE.

My question is whether there is a possibility to output the streamfunction from the pdetool or to calculate it in another way from the potentialfunction.

Any suggestion is very welcome!

Uwe

Answer 1

[Deepak Ramaswamy]

Hello Uwe,

I apologize - I had posted an answer before but that was very wrong so I deleted it :) I am not seeing any easy way to do what you want. I suggest an approach below that is a little cumbersome but maybe it is of some help. I've posted this question internally and it is possible there could be a better solution.

The idea is to solve for the potential function first, followed by the stream function problem that takes it's boundary conditions from the solution of the first problem. The first part is done is pdetool and the second in the command-line

pdetool

1. Solve Laplace problem for potential, phi, with the appropriate boundary conditions

2. Export geometry, boundary condition matrices and solution to command-line

In command-line, solve a new Laplace equation problem for psi

3. calculate the gradient for phi on the boundaries useing pdegrad and reverse and flip the signs as in

dpsidx = -dphidy and dpsidy = dphidx.

Calculate the RHS of the Neumann conditions ("g" in PDE Toolbox doc) for the new Laplace problem as in

n.grad psi = (-n1*dphidy + n2*dphidx)

4. Set Neumann conditions everywhere as in 3. except setting psi at some node to be some value (e.g. zero) to resolve the singularity. Now the PDE for psi can be solved to get the solution of the streamfunction anywhere in the geometry

Let me know if this works for you,

Deepak

[Uwe]

Hello Deepak,

Thank you very much for the reply!

Your suggestion sounds very interesting and promising. Attempting to follow your advice, after solving the Laplace equation for phi, I calculated

[dphidx, dphidy] = pdegrad(p, t, phi); % produces the gradient of psi on each triangle

dpsidx=-dphidy; dpsidy=dphidx; % to reverse and flip the signs

This is straightforward. Nevertheless in the next step specifying the Neumann conditions ( n.grad psi = (-n1*dphidy + n2*dphidx) ) I face the problem of calculating (n1, n2) that is the outward unit length normal of the boundary. This needs to be done on the triangle mesh! And the projection (-n1*dphidy+n2*dphidx) also has to work on the triangle mesh. While the concept of your solution sounds clear, I think the calculation of the Neumann boundary conditions is very difficult. May be I overlooked something … Regards Uwe

Answer 2

[Uwe]

Hello Deepak,

thank you for the useful suggestions and guidance with my problem!

Best regards

Uwe