Skip to content

Methods and Approaches

The SPECULAR project develops and integrates advanced numerical methods to enable real-time simulation of needle-tissue interactions with haptic feedback.

Overview

Our methodology combines several key techniques:

  • Model Order Reduction

Accelerate global deformation computation using reduced basis methods

  • Hybrid Formulation

Combine reduced and full models for efficient local detail capture

  • Constraint-Based Contact

Stable Lagrange multiplier formulation for contact mechanics

  • Haptic Coupling

Energy-consistent coupling for stable force feedback

Model Order Reduction

Proper Generalized Decomposition (PGD)

We use PGD to create a reduced basis that captures the essential deformation modes of the liver:

\[\mathbf{u}(\mathbf{x}, t) \approx \sum_{i=1}^{n} \alpha_i(t) \mathbf{\phi}_i(\mathbf{x})\]

Where: - \(\mathbf{\phi}_i\): Spatial basis functions (offline computed) - \(\alpha_i(t)\): Time-dependent coefficients (online computed) - \(n \ll N\): Reduced dimension vs. full FEM dimension

Benefits: - Speedup of 10-100x compared to full FEM - Preserved accuracy for global deformations - Compatible with non-linear material models

Hyper-Reduction

For non-linear materials, we employ hyper-reduction (ECSW - Empirical Cubature Hyper-Reduction):

  • Reduces evaluation of internal forces to a subset of elements
  • Maintains accuracy while reducing computational cost
  • Essential for real-time performance

Hybrid Simulation Framework

Concept

The liver is partitioned into:

Region Model Purpose
Global Reduced model Fast approximation of far-field deformation
Local Full FEM Accurate near-field around needle
Interface Coupling constraints Seamless transition

Coupling Method

The coupling ensures continuity at the interface:

  • Kinematic constraint: Displacement continuity
  • Force equilibrium: Traction continuity
  • Energy consistency: No artificial energy generation

Dynamic Update

As the needle moves, the local region updates: 1. Track needle tip position 2. Identify elements in local region 3. Update coupling constraints 4. Maintain simulation stability during transitions

Needle-Tissue Contact

Constraint Formulation

We use a Lagrange multiplier approach for contact:

\[\begin{cases} \mathbf{M} \ddot{\mathbf{u}} + \mathbf{C} \dot{\mathbf{u}} + \mathbf{f}_{int}(\mathbf{u}) = \mathbf{f}_{ext} + \mathbf{G}^T \mathbf{\lambda} \\ \mathbf{g}(\mathbf{u}) \geq 0, \quad \mathbf{\lambda} \geq 0, \quad \mathbf{\lambda}^T \mathbf{g} = 0 \end{cases}\]

Where: - \(\mathbf{g}\): Gap function (signed distance) - \(\mathbf{\lambda}\): Contact forces - \(\mathbf{G}\): Constraint gradient matrix

Implementation

Using SOFA (Simulation Open Framework Architecture): - Built-in constraint solver - Support for bilateral and unilateral constraints - Efficient matrix assembly

Sliding Contact

For needle advancement: - Contact point moves along needle shaft - Continuous update of contact constraints - Friction modeling (Coulomb model)

Haptic Coupling

Coupling Architecture

Haptic Device <--> Coupling Layer <--> Mechanical Simulation
     (1kHz)            (Stability Filter)        (Variable rate)

Energy-Consistent Coupling

Our coupling scheme ensures passivity:

\[E_{delivered} = \int_{t_0}^{t} \mathbf{f}(\tau) \cdot \mathbf{v}(\tau) d\tau \geq 0\]

This prevents energy generation that could cause instability.

Implementation Details

  • Update rate: 1kHz for force computation
  • Interpolation: Smooth force signals
  • Saturation: Force limits for safety
  • Rendering: Phantom Omni/Sensable devices

Software Implementation

SOFA Framework

All methods are implemented as plugins for SOFA:

  • ModelOrderReduction: Reduced basis computation and application
  • SoftRobots: Constraint-based modeling
  • SoftRobots.Inverse: Force control
  • Haptic: Haptic device integration

Code Structure

SOFA/
├── plugins/
│   ├── ModelOrderReduction/
│   ├── SoftRobots/
│   └── Specular/
│       ├── scene/
│       ├── python/
│       └── config/
└── applications/
    └── specular-simulator/

Validation Tests

  • Unit tests: Individual components
  • Integration tests: Full pipeline
  • Performance benchmarks: Real-time verification
  • Accuracy validation: Comparison with reference solutions

Material Models

Hyperelasticity

For liver tissue, we use a Neo-Hookean material model:

\[\Psi = \frac{\mu}{2}(I_1 - 3) - \mu \ln J + \frac{\lambda}{2}(\ln J)^2\]

Parameters identified from experimental data: - Shear modulus \(\mu\): ~1-10 kPa - Bulk modulus \(\lambda\): ~10-100 kPa

Viscoelasticity

For rate-dependent behavior:

\[\mathbf{\sigma} = \mathbf{\sigma}_{elastic} + \mathbf{\sigma}_{viscous}\]

Using Prony series or fractional derivative models.

Numerical Integration

Time Stepping

  • Implicit integration: Newmark-β or Generalized-α
  • Adaptive time step: For stability/performance trade-off
  • Sub-stepping: For haptic coupling

Solver

  • Linear solver: Conjugate Gradient, LU decomposition
  • Non-linear solver: Newton-Raphson with line search
  • Constraint solver: Projected Gauss-Seidel

Performance Optimization

GPU Acceleration

  • CUDA kernels for matrix-vector products
  • Parallel assembly of internal forces
  • Collision detection on GPU

Multi-threading

  • OpenMP for parallel element computations
  • Task-based parallelism in SOFA
  • Load balancing for hybrid models

Results

Component Performance Target
Reduced model 500Hz+ Real-time ✓
Full local model 100Hz+ Real-time ✓
Haptic coupling 1kHz Stable ✓
Visual rendering 60Hz Real-time ✓