Numerical Calculation of the Fluid Stress Tensor

I have been looking into simple ways of calculating numerically the total fluid stress tensor $\pi_{ij}$ if the flow field is known in every $x,y,z$ point. Applying the fluid stress equation (see below) isn't really straightforward since it requires an analytical velocity field. I understand this question might be "naive" but could you refer me to proper literature?

$$\pi_{ij}=-p\textbf{I}+\tau_{ij} =-p\textbf{I}+\mu(\nabla{\mathbf{u}}+\nabla{\mathbf{u}^T})$$

• Have you tried using finite differences to compute the gradients?
– Paul
Aug 8 '17 at 16:38

Where come the $p$ distribution and $\vec{v}$ field from? If you have any mesh, you can always calculate the polynomial that interpolate the variables and you can differentiate it analytically, for an arbitrary set of points. To illustrate this, $e.g$ the velocity can be approximated as $$\vec{u} = \sum{\vec{u}_i\,\phi_i(x,y,z)}$$ and therefore its gradient is $$\vec{\textrm{grad}}\,\vec{v} = \vec{\textrm{grad}}\sum{\vec{u}_i\,\phi_i(x,y,z)}=\sum{\vec{u}_i\otimes\vec{\textrm{grad}}\phi_i(x,y,z)}$$

You can create the polynomial basis such that $\phi_i(x,y,z)\in C^1(\Omega)$, $i.e$ must be continuous and derivable once (at least). Hermite polynomials are perfect for this task.

• Hmmm, let's say its the output of an experiment, e.g. PIV or some other of flow field calculation and I have all velocities for a meshed cube. Actually, it will be either FEM or FVM, but I know that this would add complication, so for now I would like to treat the output velocity and pressure field as just discrete points in space with known fluid velocity and pressure! Aug 8 '17 at 14:05
• Wouldn't it better to not try to approximate/polynomial interpolate the velocity field and then differentiate it analytically, but instead differentiate it numerically from the original $u_{i}=(x_i,y_i,z_i)$ for each point of the mesh $i$ Aug 8 '17 at 14:07
• Do you have the explicit solution for $p=p(x,y,z)$ and $v=v(x,y,z)$? You can compute easily the derivatives numerically, for, say $\phi(x)'\approx (\phi(x+\delta)-\phi(x-\delta))/(2\delta)$ with $\delta\to 0$
– HBR
Aug 8 '17 at 14:14
• That's central differences isn't it? Aug 8 '17 at 14:15
• Correct. You can set $\delta$ as small as you want if you have the analytical results.
– HBR
Aug 8 '17 at 14:16

Typically your sampled fluid velocity data would 'live' on a specific element of the mesh. For example in 3D, the velocity field 'lives' on the faces of your volume mesh. They are scalar values on faces representing the flux of fluid through that face.

This is a good explanation of how you calculate differential quantities on the discrete mesh. http://www.geometry.caltech.edu/pubs/ETKSD07.pdf

The flux on a mesh is a 2-form and the gradient of that would be a 3-form, meaning it lives on cells of the mesh. In data it would be 3 numbers stored at the cells representing difference in flux in the x direction, y dir, and z dir. The gradient operator on the mesh is literally a difference flux1-flux2. So your numbers at a cell become [flux(x+)-flux(x-), flux(y+)-flux(y-), flux(z+)-flux(z-)]. I assume you're on a regular cube mesh... If not, then everything is a little different... This is an implicit representation of gradient of u. Which is a diagonal matrix in the basis representing the orientation of the cell.

Presumably your pressure data also lives on cells of the mesh so you can simply multiply by the identity and add them to your grad(u)+grad(u)T. Your result in data would be a 3x3 matrix at each cell of the mesh.

This is an implicit representation of your stress tensor. I believe the resulting diagonal matrix will represent stresses in the orientation of the cell. The choice of discretizing velocity as fluxes forces your stress matrix to be diagonal. You can rotate the basis of the stress tensor to something else, if needed, then it won't necessarily be diagonal.

If for some reason the data you recorded for velocity doesn't 'live' on faces, that's still ok. You can do an additional step to convert from wherever your data lives to the faces. For example, if you chose to represent velocity data as (x,y,z) vectors at the cells of the mesh, you can integrate the implied velocity field across the faces, resulting in the flux representation of your velocity field. You would need an interpolation scheme to do the integration, like HBR said. But a super simple method could just be to take the difference between the component of your (x,y,z) velocity in the direction of the cell face, with the component of the (x,y,z) velocity on the other side of the cell face. This is clearly not a smooth interpolation though. There's lots of interpolation schemes.

This is done by explicit calculation. You might first reconstruct gradients using least square method, e.g. one described in paper S. Muzaferija D. Gosman - Finite-Volume CFD Procedure and Adaptive Error Control Strategy for Grids of Arbitrary Topology, JCP, 138(2), pp. 766-787, 1997. link

From there, the calculation is straight forward.