Inviscid Incompressible Flow in FreeFEM

[airwarriorg91 (Gaurav Gupta)]

Hello !
I am trying to solve the incompressible and invsicid flow in FreeFEM for a 3D geometry. I have written the code to work in parallel using PETSc. The mesh is imported from GMSH mesh file. I am not sure what is wrong with my code as the solution is not correct. I thought the dt might be a problem but changing dt doesn’t affect the result.

The weak form of the GE I used,

Boundary Conditions:

1. Dirchlet B.C. for velocity at inlet.
2. Neumaan B.C. for pressure at outlet
3. Periodic B.C. on the top, bottom, front and back walls of the domain.

varf formulation of the problem

varf GE([u1,u2,u3,p],[v1,v2,v3,q])
   =int3d(Th)(u1*v1 + u2*v2 + u3*v3
   +(un1*dx(u1) + un2*dy(u1) + un3*dz(u1))*v1*dt
   +(un1*dx(u2) + un2*dy(u2) + un3*dz(u2))*v2*dt
   +(un1*dx(u3) + un2*dy(u3) + un3*dz(u3))*v3*dt
   -(dx(v1) + dy(v2) + dz(v3))*p*dt
   -(dx(u1) + dy(u2) + dz(u3))*q*dt)
   +int3d(Th)(un1*v1 + un2*v2 + un3*v3)
   +on(1, u1=1, u2=0, u3=0)
   +int2d(Th, 2)(0*q);

The edp file for your reference. cyl_new.edp (2.5 KB)

The velocity field produced from the code. It doesn’t look correct.

Thank you for the help !

Regards,
Gaurav.

[prj]

If you have a periodic fespace, you must use the additional macro ThPeriodicity, see FreeFem-sources/examples/hpddm/diffusion-periodic-2d-PETSc.edp at master · FreeFem/FreeFem-sources · GitHub.

[airwarriorg91 (Gaurav Gupta)]

Thanks for the reply ! I will check once, is the missing macro causing the issue ?

[prj]

It is definitely one issue. I haven’t checked the rest of the code.

[cmd (Chris)]

I assume you are aware of the problems of singularities and finite-time blow-up in the incompressible Euler equations. Have you already checked your formulation with a simpler 2D/3D sequential problem to rule out such issues?

[airwarriorg91 (Gaurav Gupta)]

I was actually checking the code with the cylinder case. The problem was there in the sequential case too. No, I am not very well aware of the issues with them. Can you please share ?

[cmd (Chris)]

Well the system is hyperbolic except for the divergence-free constraint induced by the pressure. Also, the BCs in the Euler equations are not equivalent to those of the NSE (eg your no-slip condition on the cylinder does not make sense).

What is it precisely that you need from your model?

[airwarriorg91 (Gaurav Gupta)]

Okay. But, if no viscosity is present then there won’t be any “no-slip” bc on the cylinder body. There will be stagnation point in front of the cylinder but doesn’t necessarily require the velocity to be zero everywhere on the surface. So, it is free slip only, right ? (I haven’t applied any bc for velocity/pressure on the surface of the cylinder).

[cmd (Chris)]

Sorry, I only glanced at your code and misinterpreted your BC. You have not imposed no-slip BC, but you also have not imposed no-penetration BC’s… You just have a hole in your mesh.

Note that there is no natural Neumann condition in the Euler equations since there is no viscous term to integrate by parts…

[airwarriorg91 (Gaurav Gupta)]

Thanks a lot for looking into. Is it okay if put U•n =0 on all the boundaries except the inlet ?

[cmd (Chris)]

Assuming you mean all boundaries except inlet AND outlet, yes. but that is not trivial with curved boundaries. What are you trying to do? I thought you wanted periodic boundaries?

[airwarriorg91 (Gaurav Gupta)]

Yeah with curved boundary, it is any issue. I have periodic bc on the other boundaries except the cylinder boundary, inlet and outlet, does that prevent the penetration. Any help on implementing the no penetration bc for general boundaries ?

[aszaboa]

You could just solve Poisson’s equation for the complex potential; then, obtain the inviscid velocity via differentiation.

[airwarriorg91 (Gaurav Gupta)]

Stream function method doesn’t work for 3D. I have 3D application.

[cmd (Chris)]

See the notes of @fb77 here:

Boundary condition

If you want to avoid using a Lagrange multiplier (because it enlarges the matrix to invert), some other ways to penalize the normal boundary conditions are possible. In 3D you can use for example the boundary vertices. You have to define first a continuous approximate normal on the boundary (label 1 is for boundary) fespace Vh1(Th,[P2,P2,P2]); varf bdryn([w1,w2,w3],[v1,v2,v3])=int2d(Th,1)(v1*N.x+v2*N.y+v3*N.z); Vh1 [normalappx,normalappy,normalappz]; normalappx[]=bdryn(0,Vh1); Then you can ad…

[fb77 (François Bouchut)]

Hello,
The no penetration boundary conditions is one possible issue.
I think that another one is the implicit treatment of convection without viscosity. I would recommend using convect().

[airwarriorg91 (Gaurav Gupta)]

I had been using convect() earlier but it is not recommended with periodic boundary conditions. I read one discussion about it in the community. So, I switched to this method.

[fb77 (François Bouchut)]

I understand your reason, but I think that the implicit method is unstable here. You could try with convect and no penetration condition instead of periodic, to see if it works or not.

[airwarriorg91 (Gaurav Gupta)]

Sure, I will check once. Thank you for the help !

[airwarriorg91 (Gaurav Gupta)]

Hello,
Based on the recommendations, I modified my code to remove the periodic boundary condition and use convect function for the non-linear term along with no-penetration BC using the method suggested by @fb77. Thanks a lot, seems like the divergence problem is solved. But the velocity field is still patchy (See the image for reference). I don’t understand why there are patches of zero velocity. Is it due to the FEM element type. I am using P2 for velocity and P1 for pressure.

I have attached the modified code for your reference.

cyl_new.edp (3.4 KB)