# Why wall shear stress calculated from LBM directly and the one calculated based on velocity profile are so different in some cases?

First of all, I hope you accept my apologizes if my question seems off topic here. But, I asked this question in ParaView forum and after a week still I did not receive any response yet, so I'm wondering if I could find an answer hopefully here.

Background: I'm using lattice Boltzmann method to do simulation of blood flow in brain arteries. Basically, those arteries are unstructured manifolds, but we approximate them by using some tiny voxels in order to have an unstructured regular grid (I know it's a bit contradictory to have an unstructured regular grid but please bear with me). I added an example of such voxelized morphology here:

We calculate wall shear stress by using this formula from Matyka et. al., specially equation 2:

$$\vec{\tau}_{WSS} = -\mu \omega \sum_{\alpha} f_{\alpha}^{neq} \vec{c}_{\alpha} \cdot \vec{n} (\vec{c}_{\alpha} - (\vec{c}_{\alpha} \cdot \vec{n}) \vec{n})$$

Where $$\vec{\tau}_{WSS}$$ is wall shear stress vector, $$\mu$$ is dynamic viscosity, $$\tau$$ is relaxation time of LBM framework, $$f_{\alpha}^{neq} = f_{\alpha} - f_{\alpha}^{eq}$$ is non-equilibrium part of particle distribution, $$\vec{c}_{\alpha}$$ is discrete velocity in $$\alpha$$th direction, and $$\vec{n}$$ is normal vector.

The way initially during voxelization we calculate normal vector of this staircase surface is that: During voxelization, we encounter the intersection of true boundary by using CGAL intersection method, then we assign the intersected normal vector of true boundary to voxels that have intersection with that cell (triangular cell of a STL file if you want to think of it). Then, based on our distance to true boundary, we take an average in each voxel of voxelized surface based on its neighbors and their normal vectors to obtain final normal vector this staircase surface:

$$\vec{n}^{i} = \frac{\sum_{\nu = 1}^{N_{neighbors}} \frac{\vec{n}_{\nu}^{i}}{d_{\nu}}}{\sum_{\nu=1}^{N_{neighbors}}\frac{1}{d_{\nu}}}$$

Where $$N_{neighbors}$$ are the number of neighbors of the voxel in the surface, $$\vec{n}_{\nu}^{i}$$ are the initial normal vectors of voxel $$i$$th $$\nu$$th neighbor from CGAL, and finally $$d_{\nu}$$ is the distance of voxel to the true boundary in non-dimensional unit.

So, we get the wall shear stress from here and it is compared to other finite element flow solvers, which shows that error is less than 10%.

On the other hand, from macroscopic point of view, we could calculate wall shear stress from velocity profile as:

$$\vec{\tau}_{WSS} = \mathbf{T} \cdot \vec{n} - (\vec{n} \cdot \mathbf{T} \cdot \vec{n}) \vec{n}$$

Where:

$$\mathbf{T} = \mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^{T})$$

We expect that if we use above formula, we should get more or less same value as LBM calculations. Note that we already know $$\mathbf{u}$$ from LBM simulation, but we want to use GradientOfUnstructuredDataSet class from ParaView to calculate wall shear stress from velocity field. We do this for 69 different brain artery geometries. Note that we don't use exact same $$\vec{n}$$ value from voxelization part here. Instead of that, we extract boundary of voxelized mesh by taking ExtractSurface class of ParaView and we use GenerateSurfaceNormals class to create normal vector of that staircase surface.

Problem:

You can see the result of these 69 different cases in this plot, where x axis is the WSS calculated from LBM and y axis is the WSS calculated from ParaView. The thing, which is really strange is that if the above procedure were completely wrong, we should see a systematic offset of datapoints, but obviously my datapoints are all over the place and seems distributed somewhat randomly:

I took, two extreme cases where I have perfect match between LBM and ParaView and where I have complete disagreement between LBM and ParaView and their WSS distribution are shown as below.

Complete disagreement cases:

LBM:

ParaView:

Complete agreement cases:

LBM:

ParaView:

Question:

I appreciate if someone has any idea what's going here. I don't have any idea why it behaves like this. For some cases, it just shows perfect agreement, but for some other cases with exact same script and code it just produce a result completely different. Am I overlooking something hidden in the GradientOfUnstructuredDataSet? Particularly, does anyone know how GradientOfUnstructuredDataSet handles boundaries to take gradient? Still, the answer of those question does not give the answer why for some cases it works quite well and for some other ones it does a pretty bad job.

• So if I understand correctly, you are using the data set and deriving it within paraview using GradientOfUnstructuredDataSet. Do you have any idea what is the technique used for the differentiation within that method? Maybe it uses Radial Basis Function or something similar which could lead to issues close to the walls notably? Since you gave an unstructured data set, it has to derive using the little information it has. – BlaB Nov 14 '19 at 10:44
• @BlaB As far as I understand it doesn't use RBF or any other fancy differentiation techniques. It just uses simple finite difference, at least for vtkVoxel which is the cell type of my unstructured grid: github.com/Kitware/VTK/blob/… – Alone Programmer Nov 14 '19 at 14:31
• Follow up question is : What happens if you have any concave voxel? Like if you have a single row of voxel. Although your fluid velocity will be non-zero, how will it calculate the derivative in the direction where you have a single voxel? Furthermore, is it using second order differences or first order at the boundary? Because it would be problematic to use second order finite differences if you have say a layer that is two voxel wide in a given direction, because you would need two neighbhors and you would only have one. I think this may be the root of your problem – BlaB Nov 14 '19 at 16:03
• Additionally, keep in mind that in LBM, the shear stress is only approximated up to first order in space, whereas your finite difference from velocity may lead to higher spatial accuracy. Although from the looks of it, your LBM results look more realistics from a CFD point of view. – BlaB Nov 14 '19 at 16:05
• @BlaB "What happens if you have any concave voxel? Like if you have a single row of voxel. Although your fluid velocity will be non-zero, how will it calculate the derivative in the direction where you have a single voxel?" I'm somewhat surprised by your argument. Qualitatively, the cases that are off completely have more concave surface, so maybe at least it could justify why for some cases it works fine but for others it calculates nonsense. "Furthermore, is it using second order differences or first order at the boundary?" I believe based on VTK code, it uses first order at the boundaries. – Alone Programmer Nov 14 '19 at 17:56