# 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:** 20\
**Page:** 2

<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 5, 2024, 12:43pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/21 "2024-06-05T12:43:59Z")

</div>

This is quite surprising. This could be the come out of turbulence.  
Probably my computation of a smoothed normal can be better normalized, as

```auto
fespace Vh1(Th,[P2,P2,P2]);
Vh1 [normalappx,normalappy,normalappz],[w1,w2,w3];
solve bdryn([normalappx,normalappy,normalappz],[w1,w2,w3])
=int2d(Th)(w1*normalappx+w2*normalappy+w3*normalappz)
-int2d(Th)(w1*N.x+w2*N.y+w3*N.z)
+int3d(Th)(1.e-8*(w1*normalappx+w2*normalappy+w3*normalappz))
;

```

but I don’t think this will change much.

---

<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 5, 2024, 1:25pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/22 "2024-06-05T13:25:59Z")

</div>

Turbulence could be a reason but I am not able to understand how does it is getting tripped. The inlet has a constant uniform flow. Is it numerical errors ?

---

<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 5, 2024, 1:35pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/23 "2024-06-05T13:35:35Z")

</div>

Turbulence is not due to numerical error. One can describe it by saying that  
when the Reynolds number becomes large, small sctuctures appear in the solution. Their size are smaller and smaller as Re is large. The Reynolds number is Re=UL/\nu, where U is the order or magnitude of the velocity, L is the typical length of the system, and \nu is the viscosity. Here there is no formal viscosity, but you can replace it by the numerical viscosity, which is proportional to the mesh size. Thus the more you refine, the more you should see some fluctuations.

---

<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 5, 2024, 1:58pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/24 "2024-06-05T13:58:10Z")

</div>

Your code still has problems. For example, the following line (48) does nothing.

```auto
+int2d(Th, 2)(0*q)

```

Also, you have not imposed any restrictions on the lateral boundaries. Again, there is no natural condition here.

You need to work in a simpler case and build up to this. Remove the hole and get a code that works on a 2D square mesh, then add complexity.

---

<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 5, 2024, 2:09pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/25 "2024-06-05T14:09:40Z")

</div>

You could also take a larger time to see if the flow around the cylinder looks correct (but you need to visualize the direction of the velocity).

---

<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 5, 2024, 3:33pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/26 "2024-06-05T15:33:35Z")

</div>

```auto
+int2d(Th, 2)(0*q)

```

At the outlet, I wanted to implement the Neumann B.C. for pressure that’s why I used the above code. Also I have imposed the no-slip boundary condition on the top, bottom, front and back (ID=9,10,11,12) surfaces of the domain along with the curved surface of cylinder (ID = 7).

I will check with a larger time step to see if the flow around the cylinder looks correct. Otherwise, I will try with a simplier 2D domain.

---

<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 5, 2024, 3:54pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/27 "2024-06-05T15:54:10Z")

</div>

writing

```auto
+int2d(Th, 2)(0*q)

```

is equivalent to writing

```auto
+0

```

This is not a Neumann BC for pressure. There is a natural pressure release condition that arises from integrating the pressure gradient term by parts, but this is not a Neumann-type condition.

---

<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 5, 2024, 4:15pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/28 "2024-06-05T16:15:05Z")

</div>

Then, what should be outlet b.c. ? Any suggestions ?

---

<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 5, 2024, 4:19pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/29 "2024-06-05T16:19:58Z")

</div>

I think that to write no boundary term is the correct BC at outlet. It means that p=0 on this boundary, ie stress free.  
On the other boundaries except the inlet, you have put the slip condition (ie no penetration u.n=0), which is good.  
At the inlet you should put only `u1=1` i.e. u.n=1, you are not allowed to put u2=0, u3=0.

---

<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 6, 2024, 4:43am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/30 "2024-06-06T04:43:34Z")

</div>

Thanks a lot ! I will try these and let you know.

---

<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 6, 2024, 9:46am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/31 "2024-06-06T09:46:19Z")

</div>

Indeed I am not completely sure if it is possible or not to impose at the inlet u2=0, u3=0. Probably in practice it does not make any difference.

---

<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 7, 2024, 4:27pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/32 "2024-06-07T16:27:35Z")

</div>

Hello @fb77 !

Can you please confirm whether the code is correct ? Is it supposed to be 1e10 or 1e-10 ?

```auto
+int2d(Th,qft=qf1pTlump,7)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
   +int2d(Th,qft=qf1pTlump,9)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
   +int2d(Th,qft=qf1pTlump,10)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
   +int2d(Th,qft=qf1pTlump,11)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
   +int2d(Th,qft=qf1pTlump,12)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz));

```

---

<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 7, 2024, 5:01pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/33 "2024-06-07T17:01:43Z")

</div>

Yes, it is correct with 1e10.  
for the lateral boundaries which are flat, instead of this complicate penalty you can use the standard on( , ui=0) with ui the normal velocity.  
The penalty is useful only on the curved boundary.  
Also it is better to use the definition of normalapp that I gave 2 days ago.

---

<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 8, 2024, 6:05am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/34 "2024-06-08T06:05:29Z")

</div>

Thank you, I will use the newer definition.

---

<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 8, 2024, 8:29am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/35 "2024-06-08T08:29:52Z")

</div>

The freeslip boundary condition is allowing the flow to accelerate to a very high speeds even at the stagnation point. See the image for reference.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/317752c306613b5862b95677142d60c7888b9d9b.png)

Updated Problem definition. I ran the code at a higher timestep to check whether the flow goes over the cylinder or not.

```auto
varf GE([u1,u2,u3,p],[v1,v2,v3,q])
   =int3d(Th)
   (u1*v1 + u2*v2 + u3*v3
   -(dx(v1) + dy(v2) + dz(v3))*p*dt
   -(dx(u1) + dy(u2) + dz(u3))*q*dt)
   +int3d(Th)(convect([un1,un2,un3],-dt,un1)*v1
   +convect([un1,un2,un3],-dt,un2)*v2
   +convect([un1,un2,un3],-dt,un3)*v3)
   +on(1, u1=1, u2=0, u3=0)
   +int2d(Th,qft=qf1pTlump,7)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))
   +on(9, u2=0)
   +on(10, u2=0)
   +on(11, u3=0)
   +on(12, u3=0);

```

I think I will start again with a 2D problem.

---

<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 8, 2024, 1:50pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/36 "2024-06-08T13:50:05Z")

</div>

Your code looks correct. I don’t understand well your image. Anyway it is a good way to well understand the 2d first. Note however that turbulence properties in 2d and 3d maybe very different.

---

<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 10, 2024, 10:42am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/37 "2024-06-10T10:42:12Z")

</div>

The Lagrange multiplier method is another way to impose the no penetration condition on the obstacle. This can be written (for the 3d problem) as

```auto
  +int2d(Th,7)((v1*N.x+v2*N.y+v3*N.z)*pb)
  +int2d(Th,7)((u1*N.x+u2*N.y+u3*N.z)*qb)
  +int3d(Th)(1.e-8*pb*qb)

```

where `pb` is a new unknown (Lagrange multiplier) and `qb` its associated test function. They should be taken in the same space as u, i.e. P2.  
It seems to be stable with the implicit convection

---

<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 11, 2024, 5:41pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/38 "2024-06-11T17:41:08Z")

</div>

Hii Everyone !  
Update: I modified my governing equations. I did some reading and found about artificial compressibility method. I implemented it in 2D for the cylinder case and I got satisfactory results except at the stagnation point. I have attached the plot of pressure distribution along the cylinder surface varied with the artificial compressibility parameter and number of points on the surface.

 ![Variation with Beta](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/b/b27266b7971a2232c263f812632eb650e93a847c.png) ![Variation with number of points on the surface](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/9/98ffa1cfb846197ebcf0eddb6c60c0b0dd51a8e3.png)

The steady state flow looks something like this.

 ![n150](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/f/f621cb38d8f4dc325049aae09ea2b27d6c09114d.jpeg)

I will try to extend the code to 3D and check if I am getting the same results. Any idea why the pressure at the stagnation point is lower ?

---

<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 11, 2024, 10:52pm UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/39 "2024-06-11T22:52:15Z")

</div>

My guess would be boundary conditions… have you tried increasing the size of your domain? (Or decreasing the size of the cylinder)

---

<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 12, 2024, 5:22am UTC](https://community.freefem.org/t/inviscid-incompressible-flow-in-freefem/3284/40 "2024-06-12T05:22:33Z")

</div>

I have used the same boundary conditions as we had discussed earlier. Dirchlet at inlet, no b.c. on outlet and rest of the walls with slip b.c. I will try increasing the domain and check whether it helps.

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

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