Computing normal and tangent stress at boundary in contact problem with hyperelastic material

dolfinx
,
1

Hi everyone

I am solving a Nitsche-based Coulomb friction contact problem in DOLFINx for a simple triangle in contact with a flat, rigid plane (y=0). I am trying to extract and verify the normal and tangential stress components along the contact boundary.

With a linear elastic implementation, the following post-processing approach matches the analytical solutions perfectly. However, when transitioning to finite-strain hyperelasticity (using the First Piola-Kirchhoff stress tensor, PK1), the same extraction approach yields non-zero stresses on a completely free boundary.

Here is how I compute the normal and tangent stresses:

  1. Interpolate/project the stress tensor onto a ("DG", 0) function space.

  2. Loop over the boundary facets of interest to find the attached cell_id using the facet<-->cell connectivity.

  3. Extract the cell-centered stress tensor and compute the local components using the reference outward normal refr_n_hat and tangent refr_t_hat in NumPy:

Python

# Extract and compute components per cell/facet
req_stress   = stress_cells[cellid]                      # 2x2 tensor
req_tractn   = req_stress @ refr_n_hat                   # vector = tensor dot vector
req_sigma_n  = req_tractn @ refr_n_hat                   # scalar = vector dot vector
req_tractn_t = req_tractn - req_sigma_n * refr_n_hat     # vector = vector - scalar x vector
req_sigma_t  = req_tractn_t @ refr_t_hat                 # scalar = vector dot vector

I suspect that evaluating the stress via cell-centered DG0 spaces uses interior Gauss-point values. This corrupts the true boundary skin values, which should be zero for a free boundary. This is especially an issue in the hyperelastic case where the stress gradient is non-constant across a single element.

I tried using a very fine mesh and found that the same approach yields normal and tangent stress values that are ~ 5-6 orders of magnitude lower than when I use a coarse mesh.

Is there a standard or recommended way in DOLFINx to evaluate or project UFL stress expressions (PK1 or Cauchy) directly at the boundary facets/nodes without getting corrupted by interior cell values?

I would highly appreciate any kind of clarifications or support I can get on this.

Best regards

Bhanu