
[Sponsors] 
Implementation of higher order FD scheme near boundries 

LinkBack  Thread Tools  Search this Thread  Display Modes 
September 21, 2016, 04:16 
Implementation of higher order FD scheme near boundries

#1 
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
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 NS 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). HighOrder Fluxes for Conservative SkewSymmetriclike 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/00219991(92)900462 

September 21, 2016, 07:44 

#2 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,394
Rep Power: 67 
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.... 

September 22, 2016, 05:03 

#3 
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
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.. Last edited by pawank; September 22, 2016 at 07:27. 

September 22, 2016, 08:58 

#4  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,394
Rep Power: 67 
Quote:
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 

September 22, 2016, 12:08 

#5  
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
Quote:
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 5point 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 

September 22, 2016, 12:18 

#6 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,394
Rep Power: 67 
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... 

September 22, 2016, 12:28 

#7 
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
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 nonperiodic 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 

September 22, 2016, 12:32 

#8  
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
Quote:


September 22, 2016, 13:07 

#9 
Senior Member
Join Date: Jul 2009
Posts: 337
Rep Power: 17 
Have you performed a freestream 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.


September 22, 2016, 20:46 

#10 
Senior Member
Martin Hegedus
Join Date: Feb 2011
Posts: 500
Rep Power: 18 
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 nonconservativism. 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.


September 22, 2016, 20:49 

#11 
Senior Member
Martin Hegedus
Join Date: Feb 2011
Posts: 500
Rep Power: 18 
Oops, should have wrote d(rho*w)/dz instead of d(rho*u)/dz.


September 23, 2016, 01:12 

#12  
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
Quote:
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. 

September 23, 2016, 01:16 

#13  
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
Quote:


September 23, 2016, 01:20 

#14  
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
Quote:


September 23, 2016, 03:33 

#15  
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,394
Rep Power: 67 
Quote:
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. 

September 23, 2016, 08:50 

#16  
New Member
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 8 
Quote:
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. 

September 23, 2016, 08:52 

#17 
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,394
Rep Power: 67 
well, along with the BC.s implementation, I suggest to check also for the numerical stability constraints ...


Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Gradient schemes and backward second order time scheme  callumso  OpenFOAM Running, Solving & CFD  4  April 19, 2021 09:36 
Implementation of 2nd order upwind scheme  jaason  OpenFOAM Running, Solving & CFD  4  February 6, 2015 17:40 
First Order to Higher Order Blending Factor  NormalVector  FLUENT  3  April 3, 2013 15:43 
Higher order formulation question  peterciaran  Main CFD Forum  1  July 27, 2012 13:47 
High order difference scheme for KIVA3v2???  LiQiang  Main CFD Forum  1  May 9, 2005 13:50 