Introduction

Hard-magnetic soft materials, usually made by embedding hard magnetic neodymium-iron-boron (NdFeB) microparticles into soft matrix-like silicone elastomers, have attracted great attention due to their various remarkable features such as the response to remote external stimuli, fast actuation, excellent flexibility, and stretchability. [1][2][3][4][5][6][7] Promising applications include soft robotics, [8][9][10][11][12][13][14] machines and actuators, [15][16][17][18][19][20] microfluidics, [21][22][23][24][25] biomedical devices 10,13,14,[26][27][28][29][30] (e.g., endovascular neurosurgery 10 and smart catheters 29 ), and multifunctional architected materials and metastructures, 8,25,[31][32][33][34][35][36][37] just to name a few. The rapid developments of the field also call for efficient and accurate modeling and simulation platforms to rationalize the design, because these smart structures usually undergo very large and nonlinear deformation in complicated working environments, e.g., confined, enclosed spaces, combined liquid and solid terrain and brain's narrow and winding vasculature. 9,14,20,38 For example, Hu et al. 9 made untethered small-scale robots that can perform various complicated forms of locomotion, such as walking, jumping, rolling, and swimming by designing a harmonic magnetization profile in a magneto-elastic, rectangular sheet. These locomotion modes require a sophisticated control scheme to maneuver the robot as highly nonlinear deformation and multiphysics couplings are inevitable. An efficient and accurate modeling platform can provide a solid foundation for identifying the optimum control scheme and reduce the experimental trial-and-error.

Among various numerical techniques, [32][33][34][35][36][37][38][39] finite element methods (FEMs) are widely used to simulate the nonlinear and active deformation of ferromagnetic materials. 32,39 For example, Zhao et al. 32 derived a finite element simulation scheme for the hardmagnetic soft materials and implemented it into the FEM software ABAQUS through a user-defined element (UEL). This robust and powerful platform has greatly promoted the computational contributions to the design of magnetic smart structures. However, it is still difficult to include the multiphysics interactions into the framework of UEL in ABAQUS, such as fluid-structure interactions. Besides, the commercial software's parallel simulation capabilities are also constrained by the setup of the computational environments and the availability of the research licenses. Lattice models are another kind of numerical techniques for simulating elastic solids that are efficient and easy to integrate with other methods for multiphysics problems. [40][41][42][43][44][45][46][47][48] For example, lattice spring gel models have been applied to study the fluid-driven motion of microcapsules on compliant surfaces. 44,49 Recently, the authors derived irregular lattice models directly from the FEM framework, 46 and further coupled the lattice model with the lattice Boltzmann method (LBM) for complicated fluid-structure interactions (e.g., the platform of OpenFSI 50 ). Although the simplicity and powerfulness of the lattice model make it very promising, there is still a lack of a rigorous derivation to include the magnetic forces.

In this work, we present a modeling and simulation platform of hard-magnetic soft materials that combines the advantages of the recently developed FEM schemes and lattice models, named 'magttice'. We show that the magnetic actuation can be incorporated into the existing lattice model as nodal forces, which can be pre-computed if the external magnetic field is uniform. We implement the magttice model into the open-source molecular dynamics package, LAMMPS, and leverage its highly parallelized simulation capability and the versatile simulation techniques for multiphysics problems, such as fluid-structure interactions, by coupling the magttice model with the LBM. We then apply the magttice model to investigate the folded deformation of an origami plate, nonlinear deformation, and swimming of small-scale robots made by using a magnetic strip. These examples indicate that the proposed magttice model can enable more efficient mechanical modeling and simulation for the rational design of magnetically driven smart structures.

Computational model and benchmark

For hard-magnetic soft materials under moderate deformation (i.e., strain less than 30%), Zhao et al. 32 have shown that the strain energy density (U) can be expressed as the summation of elastic energy (i.e., the neo-Hookean model) and magnetic energy

where G is the shear modulus, K is the bulk modulus, I 1 = l 1 2 + l 2 2 + l 3 2 for plane strain deformation is the first invariant of the right Cauchy-Green deformation tensor, l i , i = 1, 2, 3 is the principal stretches and J is the determinant of the deformation gradient tensor F, m 0 is the vacuum permeability, B ˜r and B applied denote the residual and applied magnetic flux densities, respectively. Note that B ˜r is defined in the reference configuration. The strain energy density in eqn (1) is derived for ideal hard-magnetic soft materials, 32 whose magnetic flux density is linearly related to the applied magnetic field. For more general cases, higher-order terms of the magnetic flux density should be included in the strain energy density function. 39,51,52 For large deformation, the soft matrix will exhibit strain hardening (i.e., 30% for PDMS 53 ) and thus requires modifying the elastic strain energy density in eqn (1). Since many applications utilize large geometry nonlinear deformation with small stretching or compression, we focus on the neo-Hookean model. It should also be noted that viscoelasticity is not included in the current model and has been shown to play an important role in the deformation of magnetic soft materials. 54,55 All these nonlinearities in mechanical deformation and magnetic interactions are important and interesting research topics, which will be studied in our future work. The lattice model is derived based on energy equivalence through comparing with the finite element model. In the framework of the FEM, the deformation gradient tensor F in an irregular element can be calculated as: 56

where N a (X 1 ,X 2 ,X 3 ) is the shape function for a hexahedron element, a = 1, 2,. . .,8, and i, j, k = 1, 2, 3. The energy associated with I 1 can be written as:

By rearranging the summations in eqn (3), the energy can be obtained in the lattice model (28 lattice springs in a general case Fig. 1(a)) with the form 46

where r ab = x a À x b . Comparing above two energy expressions, we can easily calculate

The volumetric energy is calculated with an averaged volumetric strain (Fig. 1),

where V and V 0 represent the lattice's volumes at the deformed and undeformed (reference) configurations, respectively. The spring constants k ab can be computed following the same procedure of calculating thestiffness matrix in the FEM. In short, the lattice model leverages FEM techniques to rewrite the strain energy inside a hexahedron element as a summation of energies associated with lattice stretching and volumetric changes. Note that the bulk modulus K is usually not directly reported in the literature, and here we just set as K = 20 G to mimic the nearly incompressible condition of polymers and soft materials. 57,58 For the magnetic energy density, it can also be precomputed as

where

@X j dV 0 ; i;j ¼ 1; 2; 3; a ¼ 1; 2;:::;8 can be considered as general nodal forces associated with the magnetic actuation (Fig. 1(a)). It is noted that eqn (6) is valid only for spatially uniform external magnetic field because additional forces should be considered due to the gradient of the field. 59 For a given external magnetic field that is characterized by the applied magnetic flux density B applied , it is interesting to note that f m is a dead force under large deformation whose direction does not change along with the deformation. Therefore, the magnetic forces can be easily incorporated into the same framework through added nodal forces in the job script of LAMMPS simulations. This enables us to directly simulate active deformations due to magnetic forces in OpenFSI, 50 which is based on LAMMPS and already has lattice models for pure mechanics structures. For temporally changing and spatially uniform external magnetic field, we will show in the following sections (Section 3) that the updates of the magnetic forces can be directly done within the LAMMPS script without re-running the integration in eqn ( 6). This will significantly simplify the simulations of the dynamic motions of magnetic robots, such as swimming in water.

Viewing the magnetic actuation as forces can also help simplify the modeling systems. Taking a ferromagnetic beam structure as an example with B ˜r = [1, 0, 0] T and B applied = [0, 0, 1] T (Fig. 1(a) and (b)), the magnetic forces applied to the interior nodes will be cancelled out by summing the forces of the adjacent elements. This will lead to zero internal forces and non-zero nodal forces at the two ends, which is consistent with the recent observations in the modeling of the magnetic actuation of ferromagnetic beams with the elastic model. 60 We further simulate the deformed configurations of the representative ferromagnetic beam (length L = 17.5 mm, width W = 5 mm, and height H = 1 mm) using the proposed magttice model (Fig. 1(c)), where we have G = 330 kPa, K = 20 G and 1 m 0 j e B r jjB applied j as a tunable parameter. As shown in Fig. 1(d), our simulation results are in excellent agreement with FEM results reported by Zhao et al. 32 3 Applications to smart structures Since the magttice model has been validated through a benchmark of beam bending simulation, we further demonstrate its These examples include origami plates and magnetic robots, which represent structures with discretely and continuously nonuniform residual magnetic flux density, respectively.

Origami plate

Miura-ori folds have been created by encoding alternating oblique patterns of ferromagnetic domains in a simple rectangular structure (Fig. 2(a)). This is a representative example of non-uniform spatial residual magnetic flux density. We adopt the same geometries (Fig. 2(a)) and material properties in a previous work, 32 where

|B applied | = 200 mT. After discretizing the origami plate into a 3D lattice, we compute the magnetic forces at each node. We then let the structure deform under magnetic forces and relax to its final equilibrium configuration by adding damping in the simulation.

As depicted in Fig. 2(b), the magttice model can successfully capture the Miura-ori folds. It is interesting to note that forces concentrate on the domain boundaries with zero values in the interior of the domains (Fig. 2(c)). This can provide guidance for developing more efficient reduced-order models, such as plates and shells with active components, to understand the ferromagnetic origami plates' mechanical behaviors with a much larger number of unit cells.

Statically deformed configuration of the magnetic strip

We next apply the magttice model to small-scale soft-bodied robots made by using magnetic strips with continuous and nonuniform residual magnetic flux density. We take the magneto-elastic, rectangular sheet in the work of Hu et al. 9 as a representative structure and first focus on the statically deformed configurations.

As shown in Fig. 3 (the central schematic), a flat sheet is placed on a solid platform and has a harmonic magnetization profile with b R = 451. We model the same robot used in their work, with length L = 3.7 mm, width w = 1.5 mm, and thickness h = 185 mm. The material properties are also obtained from the same work, such as G = 33 kPa, K = 20 G, density r = 1.86 g cm À3 and 1 m 0 j e B r j ¼ 62 kA m À1 . In the following discussions, we keep all the geometrical and material properties the same, unless otherwise noted. To search for the static configuration after applying the external magnetic field, we gradually increase |B applied | in the simulations. After having |B applied | = 20 mT, the structure will be further relaxed for a sufficiently long time through a damped dynamic simulation until the changes of displacements become negligible. A pure repulsive interaction is assigned between the sheet and platform to describe the non-penetration conditions of the real platform. Gravity is also included in the simulation, whose important role will be discussed in detail in Section 3.3. Hu et al. 9 already showed that the robot's final configurations depended on the direction of the external magnetic field and had ''C'' and ''V'' stable configurations. To demonstrate these aspects, we run eight simulations by varying the direction y from 01 to 3151. As shown in Fig. 3, a ''V'' shape is observed when y is closer to 451 and a ''C'' shape can be obtained when y is closer to 2151. Furthermore, rigid body rotation can be noticed for y = 901, 1801, 2701, and 3151. These findings are all in excellent agreement with the analytical solutions and experiments in the previous study. 9 However, the configurations for y = 1351 and y = 3151 are more complicated, which is neither a ''V'' nor a ''C'' shape. In addition, there is no rigid body rotation and two contact points between the robot and the platform.

Rigid body motion of the magnetic strip

The case with y = 1351 or y = 3151 seems an outlier in terms of the rigid body rotation and the deformation mode. It motivates us to further examine the effects of magnetic and gravity forces applied on the robot. We notice that there is a competition between the magnetic force and gravity. Focusing on the moment at point A, the magnetic force tends to generate rigid rotation to push the strip away from the ground, while the gravity creates moments to pull the strip down to the ground. Since the rigid body rotation is also coupled with the large nonlinear elastic deformation (Fig. 4(a)), we do not pursue the analytical formula of the moment applied to point A (center of the edge (x = 0, y = 0)) at a given deformed configuration. Instead, we quantify the competition between magnetic and gravitational forces (along the negative y direction) by gradually tuning the value of density and the magnitude of the applied magnetic flux density |B applied |. When the rigid body rotation happens, only one contact point exists between the magnetic strip and platform, the tangential line of the valley of the ''V'' shape can be used as an indicator of the extent of the rigid body rotation. For the normal value of gravity on earth, no rigid body rotation is observed, even for a large magnetic magnitude as 20 mT (Fig. 4(b) and (c)). When reducing the gravity by a factor of 2, rigid body rotation begins when the magnetic magnitude is larger than 6.25 mT. And the strip reaches a plateau (B601) of the rotation angle (Fig. 4(b) and (c)). In our simulations, we can set the gravity equal to 0 (e.g., the microgravity environment) and observe almost immediate rigid body rotation (Fig. 4(b)). The theoretical rotation limit (901) is reached for magnetic magnitude as small as 5 mT (Fig. 4(c)). These examples clearly show that gravity plays an important role in determining the magnetic stripe's deformation modes, which can alter the performance of the magnetic robotic jumping, turning, and swimming. Furthermore, this also calls for more systematic simulations of the nonlinear deformation of ferromagnetic structures to design more efficient magnetic robots under different working conditions, to which the newly proposed magttice model can play an important role.

Turning of the magnetic robot

Uniform external magnetic fields with a changing direction are usually necessary to maneuver robots for various forms of locomotion, such as rolling, jumping, and swimming. To model the homogeneous external magnetic field with a changing direction, we can link the new (B 0applied ) and original (B applied ) applied magnetic flux densities through a rigid body rotation matrix

where R ij is the component of the matrix R for a given rotational angle Dy (Dy 4 0 for counter-clockwise rotation), such that

Substituting eqn ( 7) and ( 8) into the magnetic force defined by eqn (8), the new forces associated with the rotating external magnetic field can be expressed as We only need to compute the magnetic forces at each node once for a given y and directly compute the new magnetic forces after rotating the external magnetic field to y + Dy based on eqn (9). It enables us to simulate the external magnetic field's rotation by just modifying the variables in the input script of LAMMPS. We demonstrate this function by using the robot turningdeformation as an example. The robot is first deformed into a ''C'' shape by setting y = 2251 and |B applied | = 16.9 mT. Then we perform three sequential rotations along the clockwise direction with Dy = À451. The robot is rotating with the external magnetic field and maintains the ''C'' shape (Fig. 5) (see Movie 1 of ESI †). Following this deformation sequence, it is very interesting to note that a ''C'' shape can be achieved for y = 901. This is different from the ''C'' shape starting from the fixed direction of the external magnetic field (Fig. 3), indicating that the structure is bi-stable.

Swimming of magnetic robots

To demonstrate this magttice platform's capability to model the multiphysics problem, fluid-structure interactions, we investigate the swimming of the magnetic robot and compare it with the experimental study in terms of the swimming gaits of the soft robot. As shown in the experiments by Hu et al. 9 and our simulations, the soft magnetic robot's asymmetric shape can alternate between the ''C''-and ''V''-shapes under periodic B with time-varying magnitude along the principal axis. When the soft magnetic robot is immersed in water, this shapechanging mechanism can push it to swim upwards and overcome gravity, which manifests a gait similar to jellyfish swimming. 15,61,62 In the simulation, we model the same robot used in the experimental study, with length L = 3.7 mm (x-direction), width w = 1.5 mm (y-direction), and thickness h = 185 mm (z-direction). It is placed in a fluid domain with size L x = 16 mm, L y = 8 mm and L z = 16 mm with center initially at (L x /2, L y /2, L z /4). To model the swimming motion of the robot in water, a fluid-structure interaction framework is introduced. Here, we briefly review the computational method used in this framework, which is detailed in our recent work. 50 The fluid is governed by the Navier-Stokes equation (NSE)

where r f , u, and p are the fluid density, velocity, and pressure, respectively. m is the dynamic viscosity of the fluid, and F b is the body force. Instead of solving the NSE directly, we use the LBM to solve the discrete Boltzmann equation that can recover incompressible NSE through Chapman-Enskog analysis, 63 due to its high efficiency and accuracy. In our simulation scheme, the D3Q19 model is used 64 (Fig. 6(a)), where the fluid particles have possible discrete velocities in 19 directions. The explicit parameter underpinning the LBM is density distribution function f i (x,t), which is associated with the macroscopic properties of fluid as

The magnetic robot is represented by a flat sheet and accounted for by the magttice model introduced in Section 2. The mechanical properties of the magnetic robot are the same as those in the experimental study (G = 33 kPa, K = 20 G, density r = 1.86 g cm À3 and 1 m 0 j e B r j ¼ 62 kA m À1 ). The coupling between the magnetic robot and fluid is fulfilled by the immersed boundary method (IBM), which has been extensively used to study the fluid-structure interaction problems. 48,50,65,66,67 In IBM, the FSI force F f is calculated as

where b is the force coefficient, u s (t) is the velocity of the structure, and u f (t) is the fluid velocity at the position where the structure locates. u f (t) is interpolated from its surrounding fluid velocity as

where d is a smoothed approximation of the Dirac delta interpolation function. To reflect the existence of the structure in the fluid, the FSI force will be spread into the surrounding fluid mesh as

The sequence of B applied is given in Fig. 6(b), where the maximum B applied , B applied max = 17 mT and frequency f = 25 Hz. Accordingly, the robot can perform locomotion when B is oscillating along the direction a = 1051 and a = 2851. In one period, the prescribed B applied sequence can make the robot demonstrate a slow recovery stroke, in which the robot changes from the ''C''-shape to ''V''-shape (0-19.5 ms as shown in Fig. 6(c)). Following is a fast power stroke (19.5-32 ms as shown in Fig. 6(c)), which takes the robot back to the ''C''-shape. From the comparison of the configuration evolution of the magnetic robot between experiment and simulation results (cf. Fig. 6(c)), we can see that the current fluid-structure interaction platform can capture the whole swimming gaits of the robot underwater. Furthermore, the asymmetry of the ''C''-shape and ''V''-shape makes the robot experiences more propulsion in power stroke than the reduction in recovery stroke. This net propulsion will propel the robot to move upwards along the z-direction. It resembles that of an actual jellyfish swimming, 61 and it is also observed in our simulation (see Movie 2 of ESI †).