I want to get the streamline function \psi of a 2d incompressible field U=(u,v) cause the isoline of \psi are the streamlines (orthogonal to U).
I’ve seen two related topics on this forum mainly using the vorticity equation, and it seems to no working.
I’m a new Freefem user because I come from ElmerFEM. Based on the Elmer manual, the “Streamline” ElmerFEM module use the equation (46.2). I don’t understand how to implement this using Freefem cause in the weak formulation the space test function of \vec{v} is twice bigger than the space of \psi
from “Another point is that your computation of the stream function…”
The important thing is the \varepsilon.
If your u is P1, you should take the stream function in P2.
If your u is P2, you should take the stream function in P3.
I tried to implement your method on a classic thermal conduction problem (laplace and dirichlet bc), and it seems quite good, but not enough to get orthogonal streamlines with respect to the temperature.
I played with the penalty coefficient from 1e-2 to 1e-12 without any success. The P1 flux field = -grad(T) is derived from a P2 scalar field, so i used a P2 streamline scalar field on Th.
As you can see on the picture, the streamline (continuous black lines are not aligned with the vectors or not orthogonal to the temperature surface field).
Hi,
Here the zip file containing the mesh and edp file
Very strange, some streamlines are very good (on the top and on left) and others (bottom)^are not aligned with the vector flux…
Legend : fluxes (on the isostreamline meshes) are back and iso temperature are in rianbow.
I cannot see the result with Paraview, but I think that the following change can improve: integrate by parts to replace two integrals by one (result is shown below).
Indeed in the int1d there could be some internal boundaries that separate the regions. These should not be taken in the formulation.
About the int1d that includes coef there could be a similar problem.
There you can eventually replace int1d by int2d. This determines the additive constant so that \int_\Omega\psi=0 instead of \int_{\partial\Omega}\psi=0.
It looks rather nice! There is always a discrepancy because the free divergence condition on the vector field (here [fluxX,fluxY]), that is satisfied only in the weak sense. This is inherent to the discretization by finite elements of the Laplace equation (I think it is the case for any numerical method).
In your picture below it is quite surprising to have lots of discrepancies. The fact is that there are some holes in the domain, which may lead to a problem of the definition of the stream function psi. In terms of differential forms, the vector field is “closed”, but maybe not “exact”. The condition to be exact writes here \int_C \vec U\cdot\vec n dl=0,
where \vec U is the free divergence vector field, \vec n is the normal to the closed loop C, and dl is the length element along C. This must hold for any closed loop in the domain, in particular for any loop around each hole. If this is not satisfied, if follows that there exist no function \psi satisfying \mathop{\rm curl}\psi=\vec U in the domain. The variational formulation always gives a solution \psi, but in this case it does not satisfy \mathop{\rm curl}\psi=\vec U.
An explicit example of such pathology is in the annulus 1<r<2 (with r=\sqrt{x^2+y^2}) with temperature T=\log r. It satisfies \Delta T=0. We have \vec U=\nabla T. There is no stream function \psi because it should be the argument of (x,y), but we cannot define it as a continuous function in the annulus. We have \int_C \vec U\cdot\vec n dl=2\pi .The solution to the variational problem “streamlines” is \psi=0.
Yes, I do confirm the discrepancy only occurs when there are holes in the geometry otherwise it works very well (I did the same example with a ground but without holes) and the streamlines are aligned to the vector field
The problem with holes arises at the level of PDEs, it is not a problem of discrete formulation. However, a way to bypass the difficulty is to solve the streamlines problem over a modified domain, as you did (but the temperature problem can still be solved on the original mesh). We have just to cut small parts so that there remains no hole.
In the example of the annulus, one can cut the domain by considering the subdomain where the angle \theta in polar coordinates is such that |\theta|<\pi-\epsilon. Solving the streamlines formulation over this subdomain will give a correct stream function (with a jump across the cut part). The key point is to allow \psi to jump somewhere when making a loop around a hole.
I do understand the differential form is not exact due to \Omega is not a star domain, and “cutting” into small parts is not an easy task for a “general” geometry.
I just wonder myself how the cfd community manage to display the velocity streamlines when dealing with domains with obstacles (cf. picture below) ?
