Plastic hexahedral FEM for surgical simulation

Abstract

Purpose

Soft-tissue manipulations, such as collecting, stretching or tearing tissue, are a common component of surgery. When too much force is applied, these manipulations result in a residual plastic deformation that surgeons should be aware of and that should be modeled by surgical simulation.

Methods

Many tissues, vessels and organs can be modeled as offsets of curved simple shapes with primary directions, e.g., radial and axial for cylinders yield a rectangular mesh whose normal offset naturally yields a hexahedral mesh that can serve as a thick shell. Other organs are easy to embed into and deform following a hex mesh. We extend existing code for the volumetric finite element method (FEM) to model tissue plasticity as hexahedral thick shells or embedded organs. Specifically, the work extends the open source Simulation Open Framework Architecture and its newest hyperelastic deformation addition, Caribou, with focus on surgical simulation. The extension factors deformation gradients into (corotational or hyperelastic) elastic factors and plastic factors and enforces volume preservation. Limits on per-element twist, twist torque, material hardening and bounds on plasticity where elements invert avoid the need for re-meshing.

Results

Our hexahedral FEM avoids the biased outcomes of asymmetric coarsely-partitioned tetrahedral FEM. Caribou’s hyperelastic FEM is extended to hex-FEM stretching plasticity. Our high-order accurate blended-vertex deformation enables coarse hex meshes to model large plastic rotational and stretch deformations without re-meshing. We compare a vertex-blended to a cell-centered piecewise constant approach; contrast plasticity based on corotational FEM and hyperelastic FEM; and test the computation under mesh refinement. The volume is preserved also for large deformations.

Conclusion

On-the-fly generated hexahedral meshes can directly be used as finite element domains for plastic deformation based on corotational or hyperelastic elasticity. The outcome is suitable for surgical simulation.

Appendix: blending deformation gradients

Eq (19) of [6] reinterprets, for a tetrahedron with vertices $\mathbf {v}^i$, the third-order accurate ‘Phong blending’ of the center deformation $\mathbf {F}_o$ and the vertex deformations $\mathbf {F}_i$ as the ‘half-gradient’ formula: $\sum _i \beta _i \bigl (\hat{\mathbf{v}}^i + \frac{\mathbf {F}_i}{2} (\mathbf {x}-\mathbf {v}^i) \bigr )$. Here $\beta _i$, short for $\beta _i(\mathbf {x})$, are the barycentric coordinates of $\mathbf {x}$, i.e. $\mathbf {x}= \sum _i \beta _i \mathbf {v}^i$. Then $\sum _i \beta _i (\mathbf {x}-\mathbf {v}^i) = 0$ and, for a constant matrix $\mathbf {F}_o$, $\sum _i \beta _i \mathbf {F}_o (\mathbf {x}-\mathbf {v}^i)=0$ so that $\sum _i \beta _i \bigl ( \hat{\mathbf{v}}^i + \frac{\mathbf {F}_i+\mathbf {F}_o}{2} (\mathbf {x}-\mathbf {v}^i) \bigr )$ is an equally valid reformulation of the ‘half-gradient’ formula. On a hex element with 8 vertices $\mathbf {v}^{ijk}$ and $\beta _q$ the univariate barycentric coordinates of tri-linear elements $\sum _j \sum _k \sum _l \beta _j \beta _k \beta _l (\mathbf {x}-\mathbf {v}^{ijk}) = 0$ and our tri-hex analogue is $\sum _j \sum _k \sum _l \beta _j \beta _k \beta _l \bigl (\hat{\mathbf{v}}^i + \frac{\mathbf {F}_i+\mathbf {F}_o}{2} (\mathbf {x}-\mathbf {v}^{ijk}) \bigr )$.

Note that In both the corotational and the hyperelastic approach big rigid body rotations are neutralized prior to averaging $\mathbf {F}_i$ and $\mathbf {F}_0$ so that singularity of $\mathbf {F}_i+\mathbf {F}_o$ due to large rotations is ruled out. In the corotational approach, big rotations are explicitly removed from the deformation gradient. In the hyperelastic approach, the Green strain tensor (E) combines ${\mathbf {R}}^T{\mathbf {R}}={\mathbf {I}}$, i.e. rotations cancel.