This is quite surprising. This could be the come out of turbulence.
Probably my computation of a smoothed normal can be better normalized, as
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.
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 ?
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.
Your code still has problems. For example, the following line (48) does nothing.
+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.
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).
+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.
writing
+int2d(Th, 2)(0*q)
is equivalent to writing
+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.
Then, what should be outlet b.c. ? Any suggestions ?
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.
Thanks a lot ! I will try these and let you know.
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.
Hello @fb77 !
Can you please confirm whether the code is correct ? Is it supposed to be 1e10 or 1e-10 ?
+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));
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.
Thank you, I will use the newer definition.
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.
Updated Problem definition. I ran the code at a higher timestep to check whether the flow goes over the cylinder or not.
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.
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.
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
+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
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.
The steady state flow looks something like this.
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 ?
My guess would be boundary conditions… have you tried increasing the size of your domain? (Or decreasing the size of the cylinder)
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.
