The wall reflection method and wall function method
for D3Q15 Lattice Boltzmann Method
Giles Richardson (GARCFD)
17th March 2025
Required arrays
The STL normals will be stored as nrm(ntri,3) when the STL file is read in.
The STL number (1 to nstl) will be stored at each cutcell in the domain cellnrm(i,j,k).
(So that we know which STL triangle and normal are associated with each cutcell)
The cutcell number cutcell(i,j,k) will also be stored at each cell - this is just an incremental number (1:ncutc) which will later be used to access the (real) weight and (integer) normal arrays. Since we dont loop through the cutcells separately, we need a safe method to refer to each cutcell number, especially when running in parallel.
background information
The mesh used is a uniform cartesian mesh. The STL geometry is used to mark cutcells. Then a floodfill method is used to fill all the live cells. The rest of the cells remain as dead cells. For convenience this is referred to as the lcd array (live/cut/dead).
The LBM method is only applied to the live and cutcells. The dead cells are ignored. The cutcells must each have at least one live cell neighbour, otherwise they will be reassigned as a dead cell.
Whether using the bounce-back method or reflection method, it will only be applied to the cutcells, not the live or dead cells.
cutcells (grey), live cells (red), dead cells (blue)
The 15 direction vectors and weights are defined as:
The 18 ‘discrete’ normals are defined as:
The 18 discrete normals look like this, where the green points represent the normal vector endpoints on a sphere
The reflection matrix is defined as follows:
(each row corresponds to each discrete normal)
Note: Rows 1-9 are mirrored in rows 10-18
The reflection matrix looks complicated but it’s been calculated by a simple fortran code, which looks at all the incoming vectors, and reflects them depending on the surface normal.
Many surface normals were considered, and thus we are able to focus only on the surface normals which give proper reflection of the incoming vectors to the outgoing vectors.
This saves having to recalculate the reflections during the flow solver.
Calculating the nearest normals and normal weights at each cutcell
Before the solver starts we need to loop through all the cutcells and calculate the nearest 3 discrete normals and the corresponding weight for each discrete normal.
They will be stored as: weights(ncutc,3) and normals(ncutc,3)
First of all we calculate the nearest 3 discrete normals for each surface normal.
Then we store these 3 discrete normals as 3 integers (1-18) at each cutcell.
Then by calculating the distances d1,d2,d3 (from the surface normal to each discrete normal) we can assign a weight to each of the 3 discrete normals:
Calculating the distance to each discrete normals
Calculating the weights for each discrete normal
First we define a min edge length mel = 0.7654. This is roughly equivalent to the distance between closest discrete normal endpoints. Or pi/4 (chord length).
By this method we end up with 3 weights which sum to 1.0, for each cutcell.
Each weight corresponds to one of the 3 nearest discrete normals
How to implement the wall reflection method in practice
Having done the groundwork in calculating the nearest normals and the weights, the wall reflection method can then be implemented as follows:
At each cutcell we look at the nearest discrete normals, and use those to calculate the reflection vectors for each of the discrete normals eg ref1 = ref18(nn1,n), for all 15 vectors.
Then we use the weights applied to fin(i,j,k) which is then added to fout(i,j,k) for each of the reflected directions.
So to summarise, at each cutcell we use the 3 normals and 3 weights, to calculate the proportion of fin that is split 3 ways by adding to fout (for each of the the 3 reflection directions).
Coding the wall reflection method at each cutcell
summary of arrays needed
stl(ntri,3,3) STL points
nrm(ntri,3) STL normals
cutcell(ni,nj,nk) cutcell numbering system (1:ncutc)
cellnrm(ni,nj,nk) cutcell STL triangle number (1:ntri)
ref18(18,15) reflection matrix
nrm18(18,3) discrete normals
normals(ncutc,3) the 3 nearest discrete normals at each cutcell
weights(ncutc,3) the 3 weights (for each discrete normal) at each cutcell
blank slide
The wall function method
Not really a wall function - but the wall reflection method by itself will still have quite low values of velocity at the wall surface, so me may want to increase the wall velocity somewhat, like you would do by using a wall function.
By default (if you set the reach or wlfn value to zero) at each cutcell it looks for the nearest live cells, and takes the average of the fin values (for each vector).
In other words, it loops through all 27 neighbour cells (3x3x3) and looks for the live cells only. From those live cells, it takes the sum of the fin values (at each vector), and divides by the number of live neighbour cells. This gives you a new fin(n) array which can be applied at the cutcell, giving you a higher wall velocity.
This gives a robust method for increasing the cutcell wall velocity significantly.
The wall function method
But (to improve this method) instead of taking the average of the nearest live cells (to calculate fin(n) at the cutcells), we can instead “reach” out into the flowfield and grab the fin(n) values and apply those at the cutcell instead. This will give us a much higher velocity at the wall, depending on how far out (how many cells) we reach.
Admittedly - its a crude method but it works - and allows the user to decide how far (or how many cells) to reach outwards into the flowfield. By increasing the wall velocity, we can have much better control over the location and extent of separation regions.
The wall function method…
Here we assign a “reach” value (say 4 cells) to the wlfn variable. This allows us to reach out in the local surface normal direction.
The local value of fin can then be updated by using the fin values from further away from the wall, as follows:
The wall function method…
The wall function method…
When reaching out into the flowfield to look for a live cell, if the geometry is messy, we might land on a dead cell or cutcell by mistake.
So in this case, we revert to the “safe” method of using the average of the nearest live cell neighbours of the cutcell.
So there is a condition in the code which says that IF (reaching out does not give a live cell) THEN (use the average of the local live cell neighbours).