# Inviscid Incompressible Flow in FreeFEM

**URL:** <https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284>\
**Category:** General Discussion\
**Created:** [June 3, 2024, 5:28pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284 "2024-06-03T17:28:10Z")\
**Posts on this page:** 13\
**Page:** 4

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 13, 2024, 11:19am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/61 "2024-06-13T11:19:06Z")

</div>

I modified my code as

```auto
solve GE([u1,u2,p],[v1,v2,q])
   =int2d(Th)
   ([u1,u2]'*[v1,v2]/dt
   +un1*(grad(u1)'*[v1,v2])
   +un2*(grad(u2)'*[v1,v2])
   -div(v)*p
   -div(u)*q
   -Beta*p*q/dt
   )
   -int2d(Th)(
   [un1,un2]'*[v1,v2]/dt
   -Beta*pn*q/dt   
   )
   +on(1, u1=1, u2=0)
   +int1d(Th,qft=qf1pTlump,5)(1e10*(u1*normalappx+u2*normalappy)*(v1*normalappx+v2*normalappy))
   +int1d(Th,qfe=qf1pE,5)(1e10*(u1*N.x+u2*N.y)*(v1*N.x+v2*N.y))
   +on(3, u2=0)
   +on(4, u2=0);

```

In the vorticity plot, there is a bow shock like structure before the cylinder and then vorticity behind the cylinder. The explanation by @fb77 and results from the modified code present the same idea.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/6d317390951ad5ad10a775f75ff3d0ca7086104b.jpeg)

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 13, 2024, 11:23am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/62 "2024-06-13T11:23:15Z")

</div>

@fb77 Sir, can you please explain how the initial discontinuity on the left boundary is creating the vorticity wave.

Also sir, as artificial compressibility is a pseudo time. Does plotting the vorticity for smaller time means I should check the result after 5-10 iterations ?

---

<div class="post-metadata">

**Author:** ![cmd](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/cmd/32/67_2.png) [@cmd](https://community.freefem.org/u/cmd)\
**Post date:** [June 13, 2024, 12:17pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/63 "2024-06-13T12:17:21Z")

</div>

This is a canonical flow and There are many results in the literature showing numerically irritational solutions using finite element methods in primitive variables. There are still problems with your method. At r a minumum: That is not the correct representation of the nonlinear term. You may be able to diagnose these problems more quickly/ easily if you formulare this as a newton iteration rather than a pseudo time stepping scheme. Note that there are no unsteady “vorticity waves” at steady state (though there could be a steady rotational component of the solution if the BC were rotational, which they are not).

I do not wish to go back and forth about this further. I wish you the best of luck and would encourage you to read up more on this subject.

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 13, 2024, 12:19pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/64 "2024-06-13T12:19:13Z")

</div>

Thank you for all the suggestions and help ! i will surely check out the newton iteration.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [June 13, 2024, 12:23pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/65 "2024-06-13T12:23:44Z")

</div>

If you look for an incompressible potential flow u=\nabla\psi with \Delta\psi=0, you can see that it cannot satisfy simultaneously that the initial velocity is zero \psi(t=0)=0 and the inlet boundary condition \partial\psi/\partial n=-1 on the left boundary for all t\>0. Indeed if it were the case, take a smooth test function \varphi(x) which support touches arbitrarily the left boundary, but does not touch the other boundaries. Then  
\int\_{\Omega}\nabla\psi\cdot\nabla\varphi=\int\_{\partial\Omega\_{in}}(-1)\varphi.  
Since this holds for all such \varphi, and all t\>0, this prevents to have \nabla\psi(t)\to 0 as t\to 0.  
As a consequence the flow cannot be potential, therefore some vorticity is generated by the left boundary.

On the contrary if you do not put a discontinuity, for example by \partial\psi/\partial n=-t on the left boundary, then it becomes possible, thus in this case the flow remains irrotational.

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 13, 2024, 1:00pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/66 "2024-06-13T13:00:48Z")

</div>

I discussed the problem with my guide. He suggested to use convective boundary condition on the outlet instead of not specifying anything. So, I implemented the following b.c  
\frac{\partial{u}}{\partial{t}} + u \frac{\partial{u}}{\partial{x}} = 0

I updated my problem as,

```auto
solve GE([u1,u2,p],[v1,v2,q])
   =int2d(Th)
   ([u1,u2]'*[v1,v2]/dt
   +un1*(grad(u1)'*[v1,v2])
   +un2*(grad(u2)'*[v1,v2])
   -div(v)*p
   -div(u)*q
   -Beta*p*q/dt
   )
   -int2d(Th)(
   [un1,un2]'*[v1,v2]/dt
   -Beta*pn*q/dt   
   )
   +int1d(Th,2)(u1*v1/dt + un1*dx(u1)*v1 +u2*v2/dt + un1*v2*dx(u2))
   -int1d(Th,2)(un1*v1/dt + un2*v2/dt)
   +on(1, u1=1, u2=0)
   +int1d(Th,qft=qf1pTlump,5)(1e10*(u1*normalappx+u2*normalappy)*(v1*normalappx+v2*normalappy))
   +int1d(Th,qfe=qf1pE,5)(1e10*(u1*N.x+u2*N.y)*(v1*N.x+v2*N.y))
   +on(3, u2=0)
   +on(4, u2=0);

```

And it worked like charm. I think the distortion at the top and bottom points on the surface of the cylinder can be corrected by increasing the number of points. I will check but overall this is a much better solution than the previous one.

 ![Pressure Distribution](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/3592c33e0061f3724a191d7ed473c5461a2b6466.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [June 13, 2024, 2:38pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/67 "2024-06-13T14:38:51Z")

</div>

Ok it probably better drains out the vorticity coming from the initial discontinuity.

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 13, 2024, 3:21pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/68 "2024-06-13T15:21:25Z")

</div>

Yeah, seems like the case.

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [June 13, 2024, 5:11pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/69 "2024-06-13T17:11:27Z")

</div>

I was mistaken previously when I said that there was some turbulence. This would be the case if we put no slip BC on the obstacle, which is not the case here.  
About vorticity generation, I think that you could avoid it by setting no discontinuity on the inlet, as for example u1=t for 0\<t\<1, u1=1 for t\>1, whatever are the BC on the outlet.

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 14, 2024, 11:06am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/70 "2024-06-14T11:06:53Z")

</div>

Thanks a lot @cmd and @fb77 for your time and advices ! I think we can close the discussion now since the issue is resolved.

The updated code in on [GitHub](https://github.com/airwarriorg91/FreeFEM-Euler) along with the updated [report](https://github.com/airwarriorg91/FreeFEM-Euler/blob/master/Test/report/report.pdf).

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [June 19, 2024, 11:58am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/71 "2024-06-19T11:58:58Z")

</div>

An additional remark: the inlet condition u2=0 (in addition to u1=1) is incompatible with the potential flow. Thus the steady state contains some vorticity. But the condition becomes compatible when the box becomes very large, thus in this limit you recover the potential flow.

---

<div class="post-metadata">

**Author:** ![airwarriorg91](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/airwarriorg91/32/2387_2.png) [@airwarriorg91](https://community.freefem.org/u/airwarriorg91)\
**Post date:** [June 20, 2024, 5:23am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/72 "2024-06-20T05:23:49Z")

</div>

Yes sir, it makes sense. As I move the inlet close to the cylinder body the flow is not normal to the inlet it is either curved upwards or downwards. Thank You sir for pointing the fact !

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [June 20, 2024, 10:37am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/73 "2024-06-20T10:37:27Z")

</div>

Indeed there is a simple argument to see that:  
Assume by contradiction that there is a (transient or steady) potential flow u=\nabla\psi, satisfying u\cdot n=0 around the obstacle, and u=(1,0) at the inlet. Denoting u=(u\_1,u\_2), since u satisfies \mathop{\rm div}u=0 and \mathop{\rm curl}u=0, it follows that the couple (u\_1,-u\_2) satisfy the Cauchy-Riemann equations. Hence u\_1-iu\_2 is holomorphic in the domain. But on the inlet we have u\_1-iu\_2=1. By the principle of holomorphic continuation (using the Schwarz reflection principle to extend the domain of definition beyond the inlet) it follows that u\_1-iu\_2=1 in the whole domain, hence u\_1\equiv 1 and u\_2\equiv 0.  
But this contradicts the condition u\cdot n=0 around the obstacle.

Next another remark is that if you put only the condition u\cdot n=-1 at the inlet (and no condition on u\_2), you loose the uniqueness. Indeed you have a steady potentiel solution (that satisfies u\_2\not=0 on the inlet), and at least another steady flow with non-zero vorticity (that satisfies u\_2=0 on the inlet).

[Previous page](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284.md?page=3)