Bonded-particle model for magneto-elastic rods

1. Problem Statement

Hard-magnetic soft materials (hMREs) are emerging as critical platforms for untethered soft robotics, metamaterials, and biomedical devices. Simulating these systems is challenging due to the complex interplay of:

• Large Deformations: Rods and beams undergo extreme bending, twisting, and buckling.
• Non-local Interactions: Long-range magnetic dipole-dipole interactions and external field actuation.
• Contact Mechanics: Self-contact (e.g., plectoneme formation) and interaction with other bodies.
• Fluid-Structure Interaction (FSI): Operation in fluid environments (e.g., micro-swimmers, artificial cilia).

Existing approaches have limitations:

• Reduced-order rod theories (Kirchhoff/Cosserat) are efficient but struggle with complex contact and multiphysics coupling.
• Standard Bonded-Particle Models (BPM), widely used in geomechanics for 2D/3D solids, have rarely been applied to 1D structures undergoing extreme deformation, particularly with magnetic coupling.
• Discrete Differential Geometry (DDG) models handle complex shapes but often lack native support for granular contact and fluid coupling in standard simulation frameworks.

The authors aim to develop a unified, discrete-element framework that seamlessly integrates large deformations, contact, long-range magnetic interactions, and fluid dynamics within a single particle-based model.

2. Methodology

The authors developed a Bonded-Particle Model (BPM) specifically for magneto-elastic rods, implemented within the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) framework.

A. Elastic Rod Formulation (Co-rotational BPM)

• Discretization: The rod is modeled as a chain of orientable spherical particles connected by elastic bonds.
• Kinematics: Each particle possesses a position and an orientation (represented by a unit quaternion).
• Symmetric Decomposition: Unlike previous BPM formulations that use one particle as a reference frame, this model introduces a central frame ($C$) obtained by averaging the orientations of the two bonded particles via quaternion spherical linear interpolation.
• Force/Torque Calculation: Relative displacement and rotation are decomposed into four modes, each associated with a linear spring:
  1. Stretching ($K_r$): Axial force.
  2. Shearing ($K_s$): Transverse force (derived from relative displacement perpendicular to the bond).
  3. Bending ($K_b$): Torque derived from the "swing" component of rotation.
  4. Twisting ($K_t$): Torque derived from the "twist" component of rotation.
• Spring Constants: The constants are derived to match Euler-Bernoulli beam theory ($K_r = EA/L$, $K_s = 12EI/L^3$, $K_t = GJ/L$, $K_b = EI/L$).
• Symmetry: The central frame ensures the formulation is symmetric under particle exchange, eliminating rotational shear artifacts common in asymmetric formulations.

B. Magneto-Elastic Coupling

• Dipole Assignment: Each particle carries a magnetic dipole moment ($\mathbf{m}$) that rotates with the particle's orientation ($\mathbf{m} = \mathbf{R}\tilde{\mathbf{m}}$).
• Interactions:
  • External Fields: Forces and torques are calculated from the Zeeman energy ($U = -\mathbf{m} \cdot \mathbf{B}$).
  • Dipole-Dipole: Pairwise interactions are computed using standard dipole-dipole potentials, enabling long-range magnetic coupling without modifying the structural bond formulation.
• Implementation: Utilizes LAMMPS' DIPOLE package and custom bond styles (bpm/corotational).

C. Dynamics and Contact

• Integration: Equations of motion are integrated using Velocity-Verlet for translation and a Richardson update for quaternions.
• Contact: A Hertz-Mindlin granular contact model prevents self-intersection. To ensure smooth contact, particle diameters are set slightly larger than the bond length.
• Damping: Bond-level damping handles high-frequency oscillations; global Langevin damping aids convergence to quasi-static equilibrium.

D. Fluid-Structure Interaction (FSI)

• Coupling: The model is coupled with a Lattice Boltzmann Method (LBM) fluid solver using the Immersed Boundary Method (IBM).
• Mechanism: Lagrangian particles exchange momentum with the Eulerian fluid grid via a regularized delta function. Penalty forces drive particle velocities toward local fluid velocities.
• Sub-cycling: To maintain stability, multiple DEM (Discrete Element Method) steps are performed per LBM time step.

3. Key Contributions

1. Unified Framework: First implementation of a bonded-particle model for magneto-elastic rods that unifies large elastic deformations, contact, and long-range magnetic interactions in a single discrete-element code (LAMMPS).
2. Symmetric Co-rotational Formulation: Introduction of a central reference frame for bond deformation that ensures permutation invariance and stability under extreme twisting and bending.
3. Multiphysics Capability: Seamless integration of magnetic actuation and fluid dynamics (LBM-IBM) without altering the core structural formulation, allowing for the simulation of complex bio-inspired systems.
4. Scalability: Leverages LAMMPS' parallel efficiency and existing infrastructure for long-range electrostatics (adapted for magnetostatics) and fluid solvers.

4. Results and Validation

The model was validated against three benchmark problems and one multiphysics application:

• Extreme Twisting (Elastic Rods):
  • Simulated the buckling and plectoneme formation of straight and naturally curved rods under extreme twisting.
  • Result: Excellent agreement with experimental data (Lazarus et al.) and semi-analytical solutions regarding critical twist angles and deflection profiles.
• Magnetic Actuation (Non-uniform Fields):
  • Simulated magnetized beams in uniform and constant-gradient fields.
  • Result: Tip deflections matched Finite Element Method (FEM) results (Yan et al.) with high accuracy, demonstrating the model's ability to handle complex magnetic field gradients.
• Dipole-Dipole Interactions (Helical Rods):
  • Simulated helical rods with internal magnetic attraction to study mechanical hysteresis.
  • Result: Captured the hysteretic behavior (contraction vs. detachment thresholds) observed in experiments, though quantitative discrepancies were noted due to discretization limits in the point-dipole approximation.
• Fluid-Structure Interaction:
  • Oscillatory Flow: Simulated elastic filaments in oscillatory flow, matching experimental amplitude and frequency data.
  • Magnetic Cilia: Simulated fluid pumping by an array of magnetically actuated cilia with metachronal waves.
  • Result: Successfully reproduced recirculating flow patterns and demonstrated that pumping efficiency decreases with increasing wavelength, consistent with experimental observations.

5. Significance and Future Outlook

• Versatility: The particle-level approach allows for easy retrofitting of additional physics (granular, thermal, chemical) without reformulating the structural mechanics.
• Limitations:
  • Convergence: The model exhibits linear convergence (lower than the quadratic convergence of high-order FEM), requiring finer discretization for high accuracy in magnetic dipole interactions.
  • Contact vs. Resolution: The granular contact model couples particle size to bond length, limiting independent refinement of contact and magnetic resolution.
  • Energy Conservation: The current formulation uses dissipative dynamics; strict energy conservation for Hamiltonian systems requires further development of the energy potential.
• Future Work: The authors plan to extend the framework to shell-like and solid 3D structures, develop bond-based contact models to decouple resolution constraints, and address energy conservation for conservative dynamics.

In conclusion, this work provides a robust, scalable, and versatile tool for simulating complex magneto-elastic systems, bridging the gap between discrete element methods and continuum rod theories for next-generation soft robotics and metamaterials.