[ander]

Hi,

I have a question regarding how to implement boundary conditions for an unstructured cell centered fvm for compressible flow.

Say you have a no slip wall boundary, and you have one layer of ghost cells. You then want the flux contribution from the velocity to be zero. In 1D: F = [rho u, rho u + p, E]^T = [0, p, 0]^T

In the book "Computational Fluid Dynamics: Principles and Applications" by Blazek, it is proposed to set u_ghost = - u_1, (assuming u_1 is the cell nearest the wall). Up to this point I see no problem.

However, if a reconstruction procedure is used, what would now happen? The book doesn't consider this as far as I'm aware. If a reconstruction was used, The flux that would be calculated at the boundary would be F(U_L, U_R), where U_L and U_R are reconstructed values at the face. If the constant reconstruction was used, the riemann solver would return a zero velocity flux contribution, but I guess it won't for reconstructed values that doesn't necessarily obey u_R = -u_L.

Should reconstruction be turned off at boundaries?

[FMDenaro]

Originally Posted by ander

Hi,

I have a question regarding how to implement boundary conditions for an unstructured cell centered fvm for compressible flow.

Say you have a no slip wall boundary, and you have one layer of ghost cells. You then want the flux contribution from the velocity to be zero. In 1D: F = [rho u, rho u + p, E]^T = [0, p, 0]^T

In the book "Computational Fluid Dynamics: Principles and Applications" by Blazek, it is proposed to set u_ghost = - u_1, (assuming u_1 is the cell nearest the wall). Up to this point I see no problem.

However, if a reconstruction procedure is used, what would now happen? The book doesn't consider this as far as I'm aware. If a reconstruction was used, The flux that would be calculated at the boundary would be F(U_L, U_R), where U_L and U_R are reconstructed values at the face. If the constant reconstruction was used, the riemann solver would return a zero velocity flux contribution, but I guess it won't for reconstructed values that doesn't necessarily obey u_R = -u_L.

Should reconstruction be turned off at boundaries?

Indeed, the Boundary is where you prescribe the conditions not compute. However, be careful that in subsonic flows you need to let one condition from the interior.

[ander]

I am not talking about farfield conditions, I’m talking about no skip wall. My question concerns how we can ensure conservativity at the wall (no mass passing through the wall)

[FMDenaro]

Originally Posted by ander

I am not talking avout farfield conditions, I’ talking about no skip wall. My question concerns how we can ensure conservativity at the wall (no mass passing through the wall)

Since the boundary condition in a FV method are prescribed by the conditions on the fluxes of the conserved variable, you set directly the zero (or the injection condition) mass flux. Of course that requires that the face of the volume lies on the physical boundary.



In case the face of your volume is not on the boundary, I used in past to code a furhter equation for the density on the remaining part of the volume that has a face on the wall. That is congruent to the fact not necessarily you have exactly drho/dn=0 on a wall.

[ander]

But I though that the point of ghost cells is that you don’t have to give any special treatment to the boundary fluxes. Am I mistaken here? If I have do give the boundary fluxes special treatment (like setting the flux contribution from the normal velocity to zero) isn’t it easier then to drop the ghost cells entitely, and specify the boundary fluxes only based on the interior cell values and the boundary condition?

(And I’m only considering the case when the boundary is at the same location as the cell face)

[FMDenaro]

Originally Posted by ander

But I though that the point of ghost cells is that you don’t have to give any special treatment to the boundary fluxes. Am I mistaken here? If I have do give the boundary fluxes special treatment (like setting the flux contribution from the normal velocity to zero) isn’t it easier then to drop the ghost cells entitely, and specify the boundary fluxes only based on the interior cell values and the boundary condition?

(And I’m only considering the case when the boundary is at the same location as the cell face)

The necessary requirement to fulfill is to prescribe directly the fluxes if you have the face of the volume on the wall.

[ander]

But if you prescribe the fluxes directly (by directly I suppose you mean «manually» in the sense that a special flux treatment is needes at the boundary), what is then the point of ghost cells? or did I misunderstand?

[FMDenaro]

Originally Posted by ander

But if you prescribe the fluxes directly (by directly I suppose you mean «manually» in the sense that a special flux treatment is needes at the boundary), what is then the point of ghost cells? or did I misunderstand?

Could you address me exactly the page you cited from Blazek? Ghost node can be introduced for specific numerical reasons strictly related to the adopted spatial discretization.

[ander]

I am in chapter 8: Boundary Conditions. Here I'm looking at section 8.2 Solid Wall and subsection 8.2.1 Inviscid flow. I'm specifically referring to the title: "Unstructured Cell-Centred Scheme" on page 273 (in my book version).

Here it says "For the case of a triangular or tetrahedral cell,
in Ref. [5] and [6] it was suggested to employ one layer of dummy cells. The
velocity components in the dummy cells were obtained by reflecting the velocity
vectors in the boundary cells at the wall. ".

On page 275 the same concept is shown for no-slip wall. Here they take u1 = -u2 (where 1 is the ghost cell index and 2 is the nearest cell in the domain)

I would also like to point out that in section 8.1 Concept of Dummy Cells, this is written:

"The purpose of the dummy cells is to simplify the computation of the fluxes,
gradients, dissipation, etc. along the boundaries. This is achieved by the possibility to extend the stencil of the spatial discretisation scheme beyond the
boundaries. As we can see in Fig. 8.1, the same discretisation scheme can be
employed at the boundaries like inside the physical domain. Thus, we can solve
the governing equations in the same way in all “physical” grid points. ."

The sentence that I put in bold suggests to me that no special treatment of boundary fluxes should be needed, i.e that in order to assign boundary conditions, you only need to specify appropriate ghost cell values and the correct fluxes should then be calculated automatically.

However, my concern is what happens when you for instance reconstruct the variables at the faces linearly. How can you then assure that u_L = -u_R which would be necessary to get a zero velocity flux at the boundary face from the flux evaluation F(U_L, U_R).

[Eifoehn4]

Originally Posted by ander

Hi,


I have a question regarding how to implement boundary conditions for an unstructured cell centered fvm for compressible flow.

Say you have a no slip wall boundary, and you have one layer of ghost cells. You then want the flux contribution from the velocity to be zero. In 1D: F = [rho u, rho u + p, E]^T = [0, p, 0]^T

In the book "Computational Fluid Dynamics: Principles and Applications" by Blazek, it is proposed to set u_ghost = - u_1, (assuming u_1 is the cell nearest the wall). Up to this point I see no problem.

However, if a reconstruction procedure is used, what would now happen? The book doesn't consider this as far as I'm aware. If a reconstruction was used, The flux that would be calculated at the boundary would be F(U_L, U_R), where U_L and U_R are reconstructed values at the face. If the constant reconstruction was used, the riemann solver would return a zero velocity flux contribution, but I guess it won't for reconstructed values that doesn't necessarily obey u_R = -u_L.

Should reconstruction be turned off at boundaries?

FV-codes (with ghost-cells) generally apply the boundary condition twice. Once for the constant values and once for the reconstructed values. Otherwise you won't be able to perform a reconstruction at all.

The implementation is straight-forward and consistent.

If you choose to set the boundary value directly by means of flux calculation the reconstruction becomes meaningless, since the value is predescribed and fixed.

Both approaches are quite common. You may interpret the first one as weak BC and the latter one as strong BC.

A more interesting question is the wall pressure.

Regards

[ander]

I guess that makes sense. You apply the constant value first in order to be able to obtain a gradient at the domain side and then apply the reconstructed value from the domain side to the ghost cell before the flux calculation, correct?

Also, I suppose that if the ghost stencil was wide enough you wouldn't need to do it twice?

And what about the pressure?

[FMDenaro]

Originally Posted by Eifoehn4

FV-codes (with ghost-cells) generally apply the boundary condition twice. Once for the constant values and once for the reconstructed values. Otherwise you won't be able to perform a reconstruction at all.

The implementation is straight-forward and consistent.

If you choose to set the boundary value directly by means of flux calculation the reconstruction becomes meaningless, since the value is predescribed and fixed.

Both approaches are quite common. You may interpret the first one as weak BC and the latter one as strong BC.

A more interesting question is the wall pressure.

Regards

I used to evaluate wall pressure by the proper thermodynamics law wherein the density at the wall is computed by means of a evolution equation written on the volume having some faces on the wall. It worked fine on unstructured grids.

[Eifoehn4]

The exact Riemann solution of a mirror state (reflective wall BC) can be calculated analytically and is only a function of pressure. You may read Toro chapter 6.3.3.

Approximate Riemann solvers will give you a slightly different result.

[sbaffini]

My experience is that, indeed, using ghost cells for boundary conditions still requires dedicated coding and the resulting code is not simpler or more flexible (in terms of bcs you can use).

Also, none of the codes I used that relied on ghost cells were able to identically give a 0 mass flux, independently from the flow state, when using an approximate Riemann solver (altough, I can't say anything about their correctness).

I personally use direct fluxes wherever I can

[FMDenaro]

There are some hidden issues in using ghost cells. For example, Blazek addressed a linear extrapolation of the velocity but it is simple to see that such an extrapolation along the normal direction will produce the assumption that diffusive term in normal direction vanishes at second order of accuracy. Of course, for viscous flows that is not correct.

[Eifoehn4]

Originally Posted by sbaffini

My experience is that, indeed, using ghost cells for boundary conditions still requires dedicated coding and the resulting code is not simpler or more flexible (in terms of bcs you can use).

Also, none of the codes I used that relied on ghost cells were able to identically give a 0 mass flux, independently from the flow state, when using an approximate Riemann solver (altough, I can't say anything about their correctness).

I personally use direct fluxes wherever I can

Originally Posted by FMDenaro

There are some hidden issues in using ghost cells. For example, Blazek addressed a linear extrapolation of the velocity but it is simple to see that such an extrapolation along the normal direction will produce the assumption that diffusive term in normal direction vanishes at second order of accuracy. Of course, for viscous flows that is not correct.

I think conservativity is definitely not an issue. Using ghost cells will always preserve conservation. Moreover, you also do not need to implement the second approach via direct flux evaluation, since this is already implicitly available using ghost cells.

If you want to directly set the flux, simply set the both states - the inner and the ghost state - to the same value. Suprise, this will give you the physical flux, since there is no jump. Every correct implemented Riemann solver should handle this.

Moreover, a consistent implementation with regard to convective + diffusive fluxes leads directly to the use of generalized Riemann solvers.
Otherwise the approximation of both fluxes is uncoupled.

[ander]

I still have a bit of a hard time understanding how the ghost cell implementation can require no dedicated coding and be conservative at the same time for anything higher than first order. I plan on trying to implement it myself soon, so perhaps I will be back with more specific questions

[sbaffini]

I can't speak for the conservation issue but, with no doubt, a dedicated loop on the ghost cells of each boundary having a different bc is just like a loop on boundary faces to directly assign a bc.

It is just a matter of grouping and numbering. If you group and number the boundary faces together with the interior ones then, yes, it probably makes sense to use ghost cells. But if you group and number them separately, the interior face loop will never meet them, and you can loop on them separately and just do whatever you want with them.

[ander]

I would just like to clarify how the flux computation for a ghost cell unstructured code would look like. Assuming one layer of ghost cells and using Riemann solvers with linearly reconstructed face values.

Also assuming that the internal faces are grouped together and that the faces of every boundary patch are grouped together.

Here is my proposed pseudocode for evaluating the flux balance:

primitive_variables  #size = N_DOMAIN_CELLS + N_GHOST_CELLS
gradient_primitive  #size = N_DOMAIN_CELLS
limiter  #size = N_DOMAIN_CELLS
flux_balance #size = N_FACES

#step 1: assign ghost values from constant reconstruction patch wise
for boundary_patch in boundary_patches:
    BC_TYPE = boundary_patch.get_BC_TYPE() 
    for face in boundary_patch:
        i = face.i #domain index
        j = face.j #ghost index
        primitive_variables[j] = calc_ghost_value(primitive_variables[i], BC_TYPE) 

#step 2: Calculate gradient inside domain. Assuming that it used a compact stencil and that only one ghost cell is needed
#also computing the limiter.
gradient_primitive = calc_gradient(primitive_variables, other required arguments)
limiter = calc_limiter(required arguments)  

#Step 3: compute fluxes inside the domain
for face in internal_faces:
    i = face.i
    j = face.j
    primitive_reconstructed_i = #calc_reconstructed_value from gradient[i], limiter[i] and vector from center i -> face_ij 
    primitive_reconstructed_j = #calc_reconstructed_value from gradient[j], limiter[j] and vector from center j -> face_ij
    Flux = riemann_solver(primitive_reconstructed_i, primitive_reconstructed_j, face.normal)
    flux_balance[i] -= Flux
    flux_balance[j] += Flux
    
#Step 4: compute fluxes at boundaries
for boundary_patch in boundary_patches:
    BC_TYPE = boundary_patch.get_BC_TYPE() 
    for face in boundary_patch:
        i = face.i #domain index
        primitive_reconstructed_i = #calc_reconstructed_value from gradient[i], limiter[i] and vector from center i -> face_ij 
        primitive_reconstructed_j = calc_ghost_value(primitive_reconstructed_i, BC_TYPE) 
        Flux = riemann_solver(primitive_reconstructed_i, primitive_reconstructed_j, face.normal)
        flux_balance[i] -= Flux

Does the general structure here seem correct, or have I gotten something fundamental wrong?

What puzzles me a bit is that it seems to be necessary with a dedicated loop over the faces in the end. If this is indeed required, I'm starting to see why some like to evaluate the boundary flux directly.

I suppose it's not common to store the reconstructed face values (denoted primitive_reconstructed_i/j), since that would require storing a lot of additional primitive variables (two for every face). However if this was done, one could loop over all faces, reconstruct primitive variables and store them, and just do one loop over all faces to compute the flux balance.

[sbaffini]

I won't go trough the actual details of your pseudocode, but that's exactly the point.

Note that you can actually store the face value from the bc in the ghost cell and use that with a 0 gradient (pre-stored for ghosts). But your face values are still determined on a patch basis.

You can use function pointers, or pass bc functions as bc subroutine arguments but, at a certain level, some if-else statement is necessary