Hello @cmd,
I tried changing the domain size but no effect on the result. All the plots almost coincide.
You can try to impose stronger slip condition on the obstacle by setting as many conditions as degrees of freedom of velocity on the boundary. In 2d for P2 it means both lines
+int1d(Th,qfe=qf1pElump,7?)(1e10*(u1*normalappx+u2*normalappy)*(v1*normalappx+v2*normalappy))// bdry vertex
+int1d(Th,qfe=qf1pE,7?)(1e10*(u1*N.x+u2*N.y)*(v1*N.x+v2*N.y))// center of bdry edge
In 3d for P2 it means the two lines
+int2d(Th,qft=qf1pTlump,7)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))// bdry vertex
+int2d(Th,qft=qf2pT,7)(1e10*(u1*normalappx+u2*normalappy+u3*normalappz)*(v1*normalappx+v2*normalappy+v3*normalappz))// center of bdry edge
Thank You ! I will try this method.
@fb77
Sir, I tried the newer boundary condition. I didnât improve the result.
I have no further idea! What do you call âpotential flowâ? Is it an analytic solution ?
Thank you for the help everyone !
One possibility might be that I havenât imposed the condition of irrotational flow here. Although, I am comparing the solution of inviscid, incompressible and irrotational flow with that of inviscid, incompressible flow. Also, I am looking at a steady state solution during visualization, so if some vortices are being formed it might not be visible. IDK I am just speculating. i will look into it.
Yeah the potential flow refers to the analytically derived pressure distribution over a cylinder i.e. c_p = 1 - 4sin^2(\theta)
The stagnation point is behind the obstacle. I think that the pressure is modified there because of turbulent effects, that occur there because of numerical viscosity. You should see these effects if you plot the velocity vector (to see its direction). It implies also that there is some vorticity.
A similar flow is
It is a 2d Navier-Stokes flow with Re=130 against an obstacle (the little black square). The inflow is on the left. The plot is vorticity.
The picture is from Second-order entropy satisfying BGK-FVS schemes for incompressible Navier-Stokes equations
I disgree. Numerical turbulence is a spurious effect that should decay with increasing grid refinement. There also should not be numerical turbulence in 2D though there could be numerical dissipation. The unconfined inviscid flow around a cylinder is irrotational, and you should not attempt to compare the inviscid flow against results from the NavierâStokes equations due to dâAlembertâs paradox.
I believe there is indeed something wrong with your formulation. My hunch is that it is an imprecision in your BCs or a problem with your artificial compressibility implementation. But your results should match the analytical ones in the limit of a large domain.
@fb77 Yes sir, seems there is some vorticity behind the cylinder. I computed the vorticity for the velocity field (See the attached image).
I am aware about the viscous flow around a cylinder and other geometries. Actually, I am writing this code to compare the DNS solution of flow around a complex geometry with the inviscid solution of same flow.
@cmd I donât think the solution of incompressible euler equations should be necessarily irrotational. And here I am not imposing any condition for irrotationality which might be the reason for the vorticity. This implies the solution has rotational nature. This is also explains why mesh refinement or domain enlargement had any effect on the solution.
Also, I am solving the Incompressible Euler equation not the NS, hence I compared the results with the potential flow solution.
Your initial and boundary conditions are rotation-free, and the flow is conservative, therefore the flow is irrotational.
Makes sense, but then what could be triggering the vorticity generation ? You said the numerical turbulence should die down with mesh refinement and enlargement of the domain.
Find the details of the formulation in the attached pdf.
Forward Euler does not enforce incompressibility. You need to use backward euler.
Okay, I didnât knew that. Why though ? The weak formulation is using backward euler only, please check once.
It is of the form
\frac{u^{n+1} -u^{n}}{dt} = f (u^{n+1})
I think I wrote forward euler by mistake in the report.
I think that incompressibility has nothing to do with the question. Irrotationality is however related. I agree with the statement that the exact solution to the (time dependent) inviscid problem is irrotational. But the issue is that there is a numerical approximation, and its effect is similar to a small viscosity (tending to zero as the mesh size tends to zero). The result is that the solution should converge weakly to the irrotational one (as the mesh size tends to zero). But for a fixed mesh the solution well compares with the solution to the NS problem (with small viscosity, thus large Reynolds number). There are some fluctuations at very small scale. It is possible that the energy contained in these fluctuations does not tend to zero as the mesh size tends to zero. It is also possible that the time dependent solution does not tend in large time to the steady one we expect.
Ok your formulation is way off â where is the nonlinear (convection) term? This is not a linear problemâŚ
I agree with you in principle about irrotationality, but the issues displayed here are not numerical ones â there are real mistakes in his approach. Also, the large Re limit does not converge to inviscid flow: that is dâAlembertâs paradox.
Sorry, I didnât include the macros in the code there. The non-linear term is there, I used the macro con for that.
Weak Formulation:
\int_{\Omega} \left( \frac{\vec{u}^{n+1}\cdot\vec{v}}{\Delta t} + (\vec{u}^{n+1} \cdot \nabla\vec{u}^{n+1})\cdot\vec{v} - p\nabla \cdot\vec{v} - q\nabla \cdot \vec{u}^{n+1}
- \frac{\beta p^{n+1} q}{\Delta t}\right) dV - \int_{\Omega} \left(\frac{\vec{u}^{n}\cdot\vec{v}}{\Delta t} - \frac{\beta p^{n} q}{\Delta t}\right) dV= 0
macro grad(u) [dx(u), dy(u)]//
macro con(un,dt) [convect([un#1,un#2],-dt,un#1),convect([un#1,un#2],-dt,un#2)]//
macro div(u) (dx(u#1) + dy(u#2))//
solve GE([u1,u2,p],[v1,v2,q])
=int2d(Th)
([u1,u2]'*[v1,v2]/dt
-div(v)*p
-div(u)*q
-Beta*p*q/dt
)
+int2d(Th)(
-(con(un,dt)'*[v1,v2])/dt
+Beta*pn*q/dt
)
+on(1, u1=1, u2=0) //inlet
+int1d(Th,qft=qf1pTlump,5)(1e10*(u1*normalappx+u2*normalappy)*(v1*normalappx+v2*normalappy)) //cylinder wall
+int1d(Th,qfe=qf1pE,5)(1e10*(u1*N.x+u2*N.y)*(v1*N.x+v2*N.y)) //cylinder wall
+on(3, u2=0) //top wall
+on(4, u2=0); //bottom wall
Please refer to this one @cmd. I am also now convinced with the arguments of @cmd that the solution should converge to potential flow solution, as he stated the incoming flow is irrotational and the flow is conservative.
These questions are difficult, but anyway it seems to me that there exists no numerical method that could produce a solution to the initial-boundary value problem without vorticity behind the obstacle.
1 Like
No, that is not writing the nonlinear term correctly. I suggest you write it out without the convect operator to see this.
I think that the exact solution to the time dependent incompressible Euler is not irrotational. Indeed you start with vanishing initial velocity, and you impose an influx on the left boundary. It follows that there is an initial discontinuity on the left boundary. It will generate a vorticity wave that enters the domain. Then the vorticity wave propagates into the domain, and eventually goes away on the right boundary. Except that some part of vorticity gets trapped behind the obstacle because of the stagnation there.
This can be checked by plotting the vorticity for smaller time, and see if we observe the entering vorticity wave.