CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Implementation of higher order FD scheme near boundries

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 21, 2016, 04:16
Default Implementation of higher order FD scheme near boundries
  #1
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
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
pawank is offline   Reply With Quote

Old   September 21, 2016, 07:44
Default
  #2
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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....
FMDenaro is offline   Reply With Quote

Old   September 22, 2016, 05:03
Default
  #3
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
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.
pawank is offline   Reply With Quote

Old   September 22, 2016, 08:58
Default
  #4
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by pawank View Post
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
FMDenaro is offline   Reply With Quote

Old   September 22, 2016, 12:08
Default
  #5
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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
pawank is offline   Reply With Quote

Old   September 22, 2016, 12:18
Default
  #6
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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...
FMDenaro is offline   Reply With Quote

Old   September 22, 2016, 12:28
Default
  #7
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
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
pawank is offline   Reply With Quote

Old   September 22, 2016, 12:32
Default
  #8
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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...
Thank you for the suggestions. I will try this out.
pawank is offline   Reply With Quote

Old   September 22, 2016, 13:07
Default
  #9
agd
Senior Member
 
Join Date: Jul 2009
Posts: 357
Rep Power: 18
agd is on a distinguished road
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.
agd is offline   Reply With Quote

Old   September 22, 2016, 20:46
Default
  #10
Senior Member
 
Martin Hegedus
Join Date: Feb 2011
Posts: 500
Rep Power: 19
Martin Hegedus is on a distinguished road
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.
Martin Hegedus is offline   Reply With Quote

Old   September 22, 2016, 20:49
Default
  #11
Senior Member
 
Martin Hegedus
Join Date: Feb 2011
Posts: 500
Rep Power: 19
Martin Hegedus is on a distinguished road
Oops, should have wrote d(rho*w)/dz instead of d(rho*u)/dz.
Martin Hegedus is offline   Reply With Quote

Old   September 23, 2016, 01:12
Default
  #12
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
Quote:
Originally Posted by agd View Post
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.
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.
pawank is offline   Reply With Quote

Old   September 23, 2016, 01:16
Default
  #13
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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...
I tried with second order discretisations, the solver still crashes only few iterations later.
pawank is offline   Reply With Quote

Old   September 23, 2016, 01:20
Default
  #14
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
Quote:
Originally Posted by Martin Hegedus View Post
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.
Thank you.. I will think about this..
pawank is offline   Reply With Quote

Old   September 23, 2016, 03:33
Default
  #15
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by pawank View Post
I tried with second order discretisations, the solver still crashes only few iterations later.

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.
FMDenaro is offline   Reply With Quote

Old   September 23, 2016, 08:50
Default
  #16
New Member
 
Pawan Kerkar
Join Date: Sep 2016
Posts: 29
Rep Power: 9
pawank is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
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.
pawank is offline   Reply With Quote

Old   September 23, 2016, 08:52
Default
  #17
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,793
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
well, along with the BC.s implementation, I suggest to check also for the numerical stability constraints ...
FMDenaro is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
First Order to Higher Order Blending Factor NormalVector FLUENT 4 November 13, 2023 07:06
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
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


All times are GMT -4. The time now is 09:26.