WSS at boundary

[mrp089]

I'm simulating steady-state ridig-wall flow in a quarter-cylinder with constant quadratic inflow (prescribed via Spatial profile file path), constant outflow pressure, and symmetry boundary conditions (zero flow perpendicular to cut-offs).

When post-processing WSS, I get the results below. I expected WSS to be more or less uniform (no variation in circumferential direction due to symmetry and minor variation in axial direction due to developing flow). However, it seems like WSS is different at the boundaries of the surface (inlet vs. outlet, interior vs. cut-offs)

It looks like the boundaries are treated differently in the WSS calculation. @vvedula22 Before digging deeper into the WSS calculation in svFSI, do you have an idea what I would need to change to get a uniform WSS in this example? Thank you!!

[vvedula22]

@mrp089 WSS is only computed at the boundaries (faces) and not in the interior (mesh). It is set to 0 in the mesh interiors. The variations you are seeing at the face boundaries (inlet vs. outlet, cut-planes) are the side-effects of projecting WSS from elements to nodes. Even though your flow field is symmetric, since WSS is projected onto the edge shared between the wall and the symmetry faces (or the inlet/outlet faces), the absolute values of WSS on these edge nodes are underestimated. On the other hand, if you were to compute WSS as an element-based quantity, you will still see a mosaic pattern due to employing linear finite elements that will give you piecewise constant gradients. We could still provide the latter (element-based WSS) as an option.

[mrp089]

Thanks for your quick reply! I understand that WSS needs to be projected from elements to nodes. During that process, I think something is wrong in svFSI. Below, I compare WSS from svFSI (top) and "manually" in ParaView from the velocity field (bottom).

The ParaView result is uniform as expected (slight variation at the inlet due to not fully developed flow). This is how I would expect the svFSI result as well: Constant in each element gets projected to constant in each node. However, it seems that svFSI is always off in the first two layers of elements from the boundary in every direction.

@vvedula22 Do you know where BPOST is causing this?

[mrp089]

I think I know why the ParaView solution is uniform: It's by default splitting the surfaces when calculating the normals. I'm pretty sure svFSI doesn't do that, so it may calculate nodal normals that are averaged over different faces in the model.

[vvedula22]

@mrp089 I am curious how you projected WSS in Paraview but the way we project in svFSI from elements to nodes is using the lumped method where we get contributions from all the surface elements connected to the node, including adjoining faces, with appropriate weighting. In your case, because the symmetry faces and the inlet/outlet faces have negligible WSS, the edges shared with these faces get lower WSS values.

If you project from elements to nodes of the same face, then you would get a result the same as the one you obtained in Paraview. However, in svFSI, we include neighborhood faces as well.

This projection is done in BPOST and the actual projection is done here followed by getting contribution from all the processes.

[CZHU20]

@mrp089 I am also very curious about how the projection is done in Paraview. Can you elaborate on what you found? Thanks.

Lumped method works fine for P1 mesh, but can cause trouble for P2 mesh due to the near-zero weights on vertices. A better option would be a global projection, but that can be costly when factoring in the communication etc.

[mrp089]

The "trick" that ParaView uses is [Surface] Splitting (activated by default). This means that all surfaces that intersect at an angle of >30 are split in two. ParaView automatically doubles up the nodes at these intersections and calculates a separate normal (bottom). I assume that turning off Splitting (top) yields svFSI behavior, i.e. averaging the normals in each node using all connected elements.

I don't think that angle-based surface splitting is very robust in general. But we could do something similar in svFSI: We split the surfaces by faces that the user has defined in the input file. That should give a nice nodal WSS field.

[vvedula22]

Interesting. You're right that we are using face normals from a lookup table stored at nodes and that these are averaged between faces. However, even if we were to compute normals at the element level in this subroutine BPOST, we will still have this issue as we perform averaging on the final quantity of interest (e.g., WSS) at the face borders without any node duplication.

I like the idea of doing face splitting for a large separation angle, but we still have to deal with node duplication and will require changing core data structures. I am not sure if it is worth the effort. However, we could implement this as a postprocessing step using Python-based scripts. I wonder how is the WSS or another relevant quantity is processed in svPost?

[mrp089]

Yeah, messing with the discretization is probably not the way to go.

What if we export WSS (in addition) as a cell quantity as @vvedula22 originally suggested? That should circumvent any problems with normal vectors and would be a more "honest", non-interpolated way of quantifying WSS. Afterward, the user could still decide to interpolate WSS to nodes on each boundary face separately.

@ktbolt Do you know how WSS is post-processed in svPost?

[ktbolt]

@mrp089 WSS is computed in the solver, not svPost. I looked in the code and it seems that a least-squares projection is used to compute the normal to the boundary at the nodes (see solver_subroutines.f: f3lhs()) but I'm not sure what is going on, can't follow the code that well.

[vvedula22]

Thanks @ktbolt. We also do L2 or least-squares type projection to project any quantity from elements to nodes. However, it appears that in svSolver, we are using the 'special lumping technique' while in svFSI, we use the row-sum technique to do lumping of the LHS matrix.

@CZHU20 It might be worth exploring the special lumping technique's performance on higher-order elements? Please see Hughes' book pages 444-445 for more details. It is mentioned in the book that special lumping avoids negative contributions on the diagonal elements of the mass matrix typically seen in higher-order elements.

@mrp089 I agree that we should output element-based WSS in the vtu files. Will include this in the next merge.

[CZHU20]

@vvedula22 I actually implemented the special lumping techniques in one of my branches long time ago. But I remember that it didn't work that well for some reason. I will revisit my implementation and see if I can improve on it.

[mrp089]

@vvedula22 Did you have a chance to work on the element-wise calculation of WSS? I'd love to work with it if you already have something running. Thank you!!

[mrp089]

I have an update on WSS calculation. I changed BPOST to use face normals instead of Gauss-point-interpolated nodal normals. Since the geometries are split by different faces anyways, this prevents wrong normals on intersecting edges (as shown above).

Here is a comparison of different WSS formulations. All are for the same example with a steady-state flow (from left to right) through a quarter-pipe with rigid walls (no-slip).

svFSI current implementation
Velocity and normals evaluated at Gauss points. Wrong (1.2e-2 instead of 7.6e-2) at the edges due to interpolating normals over different faces.

svFSI new implementation (nodal)
Normals stay constant in each surface element, velocity evaluated at Gauss points. Note that the range [7.5e-2, 7.6e-2] is very narrow since we would expect WSS to be almost uniform over the wall.

svFSI new implementation (elemental)
As nodal implementation but averaged over all Gauss points in each element.

ParaView
Should be identical to the new implementation (nodal). Slight differences are probably due to different interpolation methods in ParaView compared to svFSI.

I'm very happy with these results so far. There is, however, one remaining issue that needs to be resolved before I can open a pull request on this.

Currently, we loop over all faces (as defined in the input file) and post-process WSS on all of them. This overwrites the results in nodes and elements that are shared between faces (and creates non-zero WSS where it doesn't make sense). Ideally, we should be post-processing WSS only on walls.

To solve this, there would be a couple of solutions:

1. Do WSS post-processing only on the first face (not robust)
2. Have the user define where they want WSS to be post-processed (complicated and easily overlooked)
3. Automatically find the correct face for WSS post-processing (no-slip in fluid simulations, interface in FSI simulations)

Solution 3 would be the most elegant, but I'm not sure if this is robust and I don't know how to implement that. @vvedula22 @ktbolt @CZHU20 Do you have any ideas?

[ktbolt]

@mrp089 You might could have another Output type for WSS where you are required to give the face names, kind of a hack though.

Output: Surface {
      Faces: lumen_wall
      WSS: t
   }

Probably the best solution is to calculate and output WSS at the element level. Then WSS could be visualized however users want. You would need to create another data structure to store results at the element level, I don't think there is anything like that in svFSI.

[vvedula22]

@mrp089 thanks for doing this. I don't think we can automatically find a wall face based on the boundary conditions as it can get complicated when we include moving walls (imposed motion). Instead, we should ask the user to provide a face on which the wall shear is computed.

In general, I think we are also interested in processing other wall quantities such as fluid pressure and traction on the wall. Currently, we have some python scripts to extract these from the vtu files but it will be a good idea to include this as part of svFSI.

As we know that all these quantities are of interest only in the walls and not on any other surface, we should ask the user to provide walls_combined.vtp file to process these results. For e.g.,

Output: Surface {
   Face file path: walls_combined.vtp
#  or we could use any face name that was already declared such as,
#  Face: < name of the face >
   Pressure: t
   WSS: t
   Traction: t
}

We can have multiple of these outputs if one is interested in these quantities on more than one face. Alternately, we can just define a list of faces within the "Output: Surface" keyword.

We should then create a set of vtp files that could be named as, "resultb_<face_name>_<time_point>.vtp", saved in the results folder.

If we are processing these quantities on the faces, we don't have to export element-based quantities as the visualization will be terrible (especially for constant strain triangles). The users then have to do additional processing to project the quantities from elements to nodes. Now that you fixed the bug in BPOST.f and are using element-based normals for WSS computation on any face, we should just export nodal quantities.

These are my ideas but open to your suggestions as well.