Implementation of higher order FD scheme near boundries
 

I have written my problem in a very detailed way. A big thank you for your time and patience to attend to my question. I appreciate this a lot.



About the solver:
The solver is FD based for 3D compressible N-S equations. For convective terms, 6th order central in split form by Ducros(1) is used, and for viscous terms, a 6th order scheme suggested by Shen(2) is used. The grid is uniform and structured. And the computational as well as the physical domain is a cuboid.

Problem:
The governing equation is dQ/dx + dE/dx + dF/dx + dG/dx = dEv/dx + dFv/dx + dGv/dx.
Considering any point in the domain, because of the schemes used, dE/dx requires a 7 point stencil(3 points on either side in x direction), while dEv/dx(which has second order derivatives) requires 4 points on either side in x direction(total of 9), plus 3 on either side in both y and z directions of each of the 9 points in x direction.

This means that the solver will need 4 ghost points in all 3 coordinate directions, which is very straight forward to handle when the boundary conditions are periodic, in fact, DNS of isotropic turbulence was done successfully. But the problem comes when the boundary conditions are not periodic.

Consider a specific case of subsonic laminar boundary layer. There are inflow, wall and outflow boundaries. Z direction is dropped in this case. Out of the 4 ghost cells outer one is taken as boundary where characteristic based boundary conditions(3) are applied. The remaining ghost cells are now part of the computational domain. On these points, wherever possible, the original scheme is used, and on the remaining points simple 4th order differences(which may be obtained from Taylor series) are used. Central differences are used wherever possible, and points where central is not possible, i.e. points just next to the boundary, forward or backward differences are used taking one point on the boundary and remaining points on the other side. For viscous terms, first the stresses and heat conduction terms are computed, then their derivatives are computed.

This did not work.

My question is, if do everything correctly, should this work? Are there any inherent problems with the schemes used, especially the use of simple finite differences, and in the way of calculating viscous fluxes?




1 Ducros, F., Laporte, F., Soulères, T., Guinot, V., Moinat, P., & Caruelle, B. (2000). High-Order Fluxes for Conservative Skew-Symmetric-like Schemes in Structured Meshes: Application to Compressible Flows. Journal of Computational Physics, 161(1), 114–139. http://doi.org/10.1006/jcph.2000.6492
2 Shen, Y., & Zha, G. (2010). Large eddy simulation using a new set of sixth order schemes for compressible viscous terms. Journal of Computational Physics, 229(22), 8296–8312. http://doi.org/10.1016/j.jcp.2010.07.017
3 Poinsot, T. J., & Lelef, S. K. (1992). Boundary conditions for direct simulations of compressible viscous flows. Journal of Computational Physics, 101(1), 104–129. http://doi.org/10.1016/0021-9991(92)90046-2

It would help if you give details about your error in the solution...

If you set the proper boundary condition for compressible (subsonic/supersonic) flow, if you use proper backword/forward discretizations near the boundaries, it should work.
Not that the stability constraint can be affected by the discretization near the boundary.
Then, I am not able to thing nothing else than checking for common programming errors....

Thank you FMDenaro.

To be precise, the pressure becomes negative at the outflow boundary at a point just next to the corner point between(intersection of) the outflow boundary and the wall boundary. And from there the problem starts.

And does the fact that the problem first appears at the boundary indicate that I have made an error in implementing boundary conditions? Or can it be something wrong with implementation forward/backward discretisation also? I have been stuck with this for 5 weeks now:( so any kind of help will be appreciated...

Thanks in advance FMDenaro and anyone else who helps.. :)

I cannot give one answer but only suggest for some checks..

1) how do you set the boundary in outflow? do you expect supersonic or subsonic outflow?
2) how do you discretize the derivatives in the points before the outflow? The large stencil implies to specify some numerical BC also for points in the interior close to the outflow

1a Outflow boundary is set using characteristic based BCs(http://doi.org/10.1016/j.jcp.2010.07.017) wherein governing equations are expressed in terms of characteristic wave amplitudes. The wave amplitudes of incoming waves are calculated based on local, one dimensional inviscid(LODI) relations. The wave amplitudes of outgoing waves are calculated using one-sided descretisations using points from the interior. Then the governing equations are solved to get values of primitive variables on the boundary. This method avoids extrapolation and numerical BCs are set using LODI.
1b The inflow is mach 0.5 uniform flow, and flat plate on the bottom so I expect the outflow to be subsonic.

2 For the two points next to the boundary, I use a fourth order descretisation which means a 5-point stencil. For the point just next to boundary, one point is taken on the boundary and 3 on the other side. And for the other point, i use central difference with two points either side(one of the points being on the boundary)

Do I need any additional BCs?

Also can you please point me to a good source of reading which talks about boundary conditions relating to my problem, i.e. a subsonic and supersonic flat plate boundary layer(inflow, outflow and wall boundary) with high order scheme such that there are not enough points available for the points near boundary? May be different types of boundary conditions, i.e. not necessarily characteristic based BC.

Thanks a lot again:)

To check definitively if the high order discretization is the source of the problem I suggest to run your code first with a standard second order accurate discretization. I understand that you have to change some lines of your code but I would spend time in doing this debug and get this answer...

Further possible problems are the numerical instability of your scheme for the computational parameter you are setting, wrong implementation of the BC.s, simple bug in the code...

The post just before this one was a reply to FMDenaro, but for anyone reading this, can you please point me to a good reading source talking about non-periodic boundary conditions for a higher order FD schemes which has a problem that there would not be enough points near boundary which is required by a higher order scheme. I want to simulate a subsonic and supersonic boundary layer for now, so I will have to implement inflow, outflow and wall BCs.

Thanks in advance:)

Thank you for the suggestions. I will try this out.

Have you performed a free-stream check of your solver? If not, then simply set your boundaries to Dirichlet, initialize your solver to freestream (or uniform values) and see if the solver maintains that uniform state. With gradient boundary conditions, this can also be a useful check if you are using ghost cells, since there should be no gradient if the flow field is completely uniform. If you have done that, then I would echo FMDenaro's suggestion.

I gave this a quick read so I'm not sure if I understand the issue, but your problem may be on the surface boundary since you are implementing a no slip condition without extrapolating to the surface. The problem is that, since the surface has zero velocity (i.e. zero flux), density and energy can accumulate there, if you are implementing the equations in conservative form. If this sounds like it could be the problem, I would suggest throwing in a little non-conservativism. For example, d(rho*u)/dz = 0.99*d(rho*u)/dz + 0.01*(u*drho/dz + rho(du/dz)). U is zero so that can be dropped, but du/dz is not. Thus rho and rho*e0 feed back into the system. That was the problem I had at the stagnation points when implementing a higher order inviscid case.

Oops, should have wrote d(rho*w)/dz instead of d(rho*u)/dz. :p

Firstly, thanks for the suggestion.

My solver is actually a hybrid one with 6th order central in Ducros split form and 5th order WENO based on lax Friedrich flux splitting. When I switch on WENO, the free stream check was successful, but not for Ducros... I'll try to figure out why.

I tried with second order discretisations, the solver still crashes only few iterations later.

Thank you.. I will think about this..

Well, somehow this is a good point to work on, it is useless to continue searching for an error in the high order implementation.
Now focus on possible bugs in the implementation of the BCs.
Honestly if your run diverges in few time steps I would think rather in a bug than in a numerical stability problem, but however how do you set the time step? Are you sure to fulfill the stability constraint?

It would be useful if you do a single time step iteration and then check for all the variables. You could also post the figures to help us understanding.

Thanks.

Sorry if I may have communicated incorrectly. What I meant was, other things being same, if the code with 4th order implementation diverges at 61st iteration then the code with 2nd order diverges at 81st iteration.

well, along with the BC.s implementation, I suggest to check also for the numerical stability constraints ...