I am trying to solve the compressible, turbulent Navier–Stokes equations for a simple test case: flow around a sphere. To exploit symmetry I model half the sphere in cylindrical coordinates. I use FreeFem++ for the variational formulation and PETSc as the nonlinear solver backend.
To my eye the boundary conditions seem to be applied correctly in the code, but the solution field looks unexpected (see attached image) : we observe density drop near the inlet. I’m not sure whether the issue is:
an implementation bug (wrong BC application)
a modelling error (BCs inappropriate for compressible/turbulent model)
a numerical issue (axis singularity, initial guess)
a mesh issue
Could someone explain to me why I have this density drop at the inlet ?
For information, the residual operator is given by :
Dear Ilyan,
I give you my opinion (the following interpretation is not necessarily correct).
I think that your boundary conditions on the inlet may be the problem.
You have positive viscosity thus you can impose the three components of the velocity. I guess you have thermal conductivity thus you can impose temperature.
For density, it does not solve an elliptic problem, but you can impose it if the normal velocity ux is positive.
This is probably what you have (you set ux=u). But if you imagine to solve a time dependent problem, you need to wait for a sufficiently long time for the boundary value of the density to propagate in the whole domain (T=L/u).
Here you don’t solve a time-dependent problem, but you solve the steady flow by Newton iteration. Since the equation on the density is not elliptic, it propagates slowly, and thus maybe you need a huge number of iterations to obtain the steady flow. I can imagine that each iteration propagates the density further of only a few cells (to be confirmed).
Another possibility is the following. I guess that the Reynolds number is large. Thus it is significant to look at the status of the inflow boundary conditions in the inviscid regime. Then if you are in the subsonic regime (Mach number less than 1), it is not possible to impose both rho and ux on the boundary, because one wave is going out of the box. You can only impose one quantity, either rho or ux, or the mass flux rho*ux (which is probably the most physically relevant). When imposing both rho and ux with viscosity, this leads to a boundary layer. This may be what you observe.
Hi Ilyan,
You have too many boundary conditions and are overconstraining the flow at the inlet here, leading to an inconsistency. I would suggest eliminating any BC on the density, as that should be enforced algebraically through the equation of state.
Note: You could also eliminate the temperature BC and keep the density BC. However, typically, you would want to keep a temperature BC since you often use temperature gradient BCs on, e.g. adiabatic walls or outflow conditions.
Thank you for your replies! I’d like to set the total pressure and total temperature at the inlet. Since I also know the inlet velocity, it seems to me that by specifying both the temperature and the density (assuming a perfect gas), the state is fully defined. That’s why I don’t understand why this would over-constrain the flow. To me, the equation of state serves to determine the pressure once the density and temperature are known. Could you clarify this for me?
Because the total pressure has already been constrained by your free-outflow condition on the outlet. You cannot constrain the pressure at both the inlet and the outlet consistently with Dirichlet velocity BC at the inlet.
To simplify, consider a 1d problem where the temperature is eliminated (it has its own equation with boundary conditions).
For a steady inviscid compressible flow \partial_x(\rho u)=0, \partial_x(\rho u^2+p(\rho))=0,
where p(\rho) is a pressure law.
Then you see that \rho u=c_1 a constant, and \rho u^2+p(\rho)=c_2 another constant, so that \rho and u are both constant.
If you impose both values of \rho and u at the inlet, then c_1 and c_2 are determined, and you cannot impose further BC on the right (or conversely if you impose 1 BC on the right, you can impose only 1 BC on the left).
With viscosity \partial_x(\rho u)=0, \partial_x(\rho u^2+p(\rho)-\nu\partial_x u)=0.
Then \rho u=c_1, \rho u^2+p(\rho)-\nu\partial_x u=c_2.
In this case you need a third constant c_3 (for example u(0)=c_3) to determine the solution u to this first-order differential equation.
Next when \nu is small, two cases can occur:
If the constants c_1,c_2,c_3 are so that the solution is compatible to the inviscid solution (1 relation between the 3 constants) then as \nu\to 0 the viscous solution tends to the inviscid solution.
If the relation is not satisfied, then as \nu\to 0 the viscous solution tends to an inviscid solution but you will not get the three BC at the limit (at least one BC is lost).
In this case there is a boundary layer: the solution u_\nu varies rapidly close to the boundary.
