*A. SMA model*

The Muller-Achenbach model is used to describe the SMA material in this study [34]. Utilizing concepts rooted in statistical thermodynamics, this model clarifies the kinetics of two separate fractions of the martensite phase. Moreover, the SMA model employs the principles of thermally activated mechanisms to analyze these processes. Drawing inspiration from crystallographic observations, the proposed model conceptualizes a mesoscopic lattice layer as the fundamental building block. The composition of the SMA material consists of austenite A alongside two distinct martensitic variants, labeled as M+ and M−, and the proportions of these phases are represented as *x*A, *x**+* and *x**−*, respectively. These quantities adhere to the subsequent relationship,

$${x}_{\text{A}}+{x}_{-}+{x}_{+}=1 \left(1\right)$$

In the domain of SMA, the expression denoting the macroscopic strain can be formulated as a summation of weighted strains pertaining to the distinct phases,

$$\epsilon ={x}_{\text{A}}\stackrel{-}{{\epsilon }_{\text{A}}}+{x}_{+}\stackrel{-}{{\epsilon }_{+}}+{x}_{-}\stackrel{-}{{\epsilon }_{-}} \left(2\right)$$

where \(\stackrel{-}{{\epsilon }_{\text{A}}}\), \(\stackrel{-}{{\epsilon }_{+}}\), and \(\stackrel{-}{{\epsilon }_{-}}\) are the average strains in each phase, which are determined by Boltzmann statistics [35]. Presuming the Gibbs free energy density g to express as a multi-parabolic function procession upon the strain acting as the pertinent order parameter, the mean strains can be ascertained by assessing the abscissa values that correspond to the minima of the corresponding energy wells.,

$$\stackrel{-}{{\epsilon }_{\text{A}}}=\frac{\sigma }{{E}_{\text{A}}} \stackrel{-}{{\epsilon }_{+}}=\frac{\sigma }{{E}_{\text{M}}}+{\epsilon }_{T} \stackrel{-}{{\epsilon }_{-}}=\frac{\sigma }{{E}_{\text{M}}}-{\epsilon }_{T} \left(3\right)$$

The explicit stress-strain relationship can be obtained by inverting the combined effect of Eq. (2) and Eq. (3),

$$\sigma \left(\epsilon \right)=\frac{{E}_{\text{M}}\left(\epsilon +{x}_{-}{\epsilon }_{T}-{x}_{+}{\epsilon }_{T}\right)}{{x}_{-}+{x}_{+}+\frac{{E}_{\text{M}}}{{E}_{\text{A}}}{x}_{\text{A}}} \left(4\right)$$

Moreover, Eq. (4) will be utilized in a conventional finite element framework that is driven by displacement. The dynamics of the phase fraction \({x}_{\alpha }\), \(\alpha =⟨\text{A},+,-⟩\) is governed by the rate laws,

$$\frac{\partial {x}_{+}}{\partial t}=-{x}_{+}{p}^{+\text{A}}+{x}_{\text{A}}{p}^{\text{A}+}$$

$$\frac{\partial {x}_{-}}{\partial t}=-{x}_{-}{p}^{-\text{A}}+{x}_{\text{A}}{p}^{\text{A}-} \left(5\right)$$

where \({p}^{\alpha \beta }\) represents the transition probability from phase \(\alpha\) to phase *β*. The transition probabilities are obtained through statistical thermodynamic analysis, influenced by the combination of material stress \(\sigma\) and temperature *T*. These factors determine the probability of an optimal layer surpassing an energy screen within a non-convex free-energy topography. The progression of *T* within the SMA material is elucidated via the application of the internal energy balance principle.

$$\rho \xi \frac{\partial T}{\partial t}-k\frac{{\partial }^{2}T}{\partial {{x}_{i}}^{2}}=h\left(\frac{\partial {x}_{+}}{\partial t}+\frac{\partial {x}_{-}}{\partial t}\right)-\alpha {S}_{v}\left(T-{T}_{0}\right)+j\left(t\right) \left(6\right)$$

where \(\rho\) is the volumetric mass density of the material, \(\xi\) signifies the specific heat, \(k\) represents the thermal conductivity, \(h\) stands for the latent heat of transformations, \({S}_{v}\) denotes the ratio between external surface area and volume of the wire, and \(j\left(t\right)\) represents the volumetric density of Joule heating [36]. This study primary focus on the comprehensive behavior of the prototype. Therefore, it is reasonable to overlook the localized phenomena occurring within SMA and adopt the presumption of a uniformly distributed temperature along the wire. Consequently, the equations (5)-(6) is able to restated as decentralized one-dimensional equations within the FEA model. Furthermore, the transition probabilities are contingent upon the transformation stresses exhibited by the austenitic and martensitic phases, which are articulated as,

$${\sigma }_{\text{A}}={\sigma }_{\text{L}}+\frac{d\sigma }{dT}\left(T-{T}_{\text{L}}\right) \left(7\right)$$

$${\sigma }_{\text{M}}={\sigma }_{\text{A}}-{\Delta }\sigma \left(8\right)$$

where \({\sigma }_{\text{L}}\) represents the transition stress occurring during the transformation from the martensite phase to the austenite phase at the specified reference temperature \({T}_{\text{L}}\), \(d\sigma /dT\)denotes the rate of change of \({\sigma }_{\text{L}}\) concerning changes in temperature *T* within the context of gradient analysis, while hysteresis width is denoted as \({\Delta }\sigma\) and is defined as \({\sigma }_{\text{A}}{\sigma }_{\text{M}}\).

*B. Implementation of FEA method in SMA model*

This study utilizes the commercial platform COMSOL Multiphysics to implement the FE model to capture the kinetic information of the actuator. The incorporation of SMA mechanical dynamics into COMSOL simulations typically necessitates the utilization of the momentum conservation [25]. Nonetheless, the existing formulation is not appropriate for incorporating the mechanical linkage between SMA springs and the connection structure. Consequently, in this study, the truss element is utilized to incorporate constitutive SMA models. Employing a truss, designed to accommodate solely axial loads, effectively captures the essence of the one-dimensional SMA model. Moreover, utilizing the truss element facilitates the establishment of a simpler connection between SMA and external units. The truss stress constitutive model is formulated as follows,

$${\sigma }_{t}={\sigma }_{ti}+E\left(\frac{du}{dx}-{\epsilon }_{in}\right) \left(9\right)$$

where \({\sigma }_{ti}\) is the initial truss stress, which defines as zero, *E* represents the truss Young’s modulus, \(du/dx\) donates the truss strain, and \({\epsilon }_{in}\) expresses the truss strain. To establish a connection between the truss units and the SMA model, it is necessary to equate the truss stress \({\sigma }_{t}\) to the SMA stress \(\sigma\). This congruence can be realized by defining the following variables,

$$E=\frac{{E}_{\text{A}}{E}_{\text{M}}}{{E}_{\text{A}}\left({x}_{+}+{x}_{-}\right)+{E}_{\text{M}}\left(1-{x}_{+}-{x}_{-}\right)} \left(10\right)$$

$$\frac{du}{dx}=\epsilon \left(11\right)$$

$${\epsilon }_{in}={\epsilon }_{T}\left({x}_{+}-{x}_{-}\right) \left(12\right)$$

In addition, in the context of the SMA model, the initial state of the spring is considered to be fully martensitic, and the initial temperature is defined as the ambient temperature.

*C. Coupling strategy between FEA and MBD*

The Multibody Dynamics (MBD) approach is employed to characterize the complete motion information of the gripper, facilitating the simulation of the problem involving substantial deformations and the loading process. This module enables the simulation of displacements and rotations in rigid bodies, along with configurations involving interconnected rigid bodies through articulating joints. Concurrently, it is possible to incorporate user-defined differential equations as constitutive relationships tailored to particular elements. Specifically, the coupling between MBD model and FEA model is established to describe the information transition between the rigid units and SMA spring. The utilization of a coupling strategy facilitates the incorporation of intrinsic attributes from MBD module for the purpose of constraining the movement of the spring within the phalanxes. The finer tetrahedral elements are used in rigid structures and the standard hexahedral elements are still used for SMA springs. In addition, all the nodes of the tetrahedral mesh are assigned to be the SMA truss. The truss element is interconnected with the elastic rigid element and can be represented as follows,

$$\left({u}_{t},{v}_{t},{w}_{t}\right)=\left({u}_{r},{v}_{r},{w}_{r}\right) \left(13\right)$$

which means the degrees of freedom of the elements and solid elements are same. Furthermore, the embedded characteristic of the MBD module additionally facilitates the establishment of a rigid linkage between the two components within the context of any motion frame.

In addition, the coupling strategy is shown in Fig. 5. At each time step of FEA-MBD coupling algorithm, the stress and strain are calculated in FEA solver while the displacement is calculated by MBD. Notably, in the coupling process, the nodes associated with rigid connection are loaded with stress information from the FEA solver for displacement calculation. Besides, the presence of different rigid connections restrict the structure movement to some extent, thus MBD also depicts the motion characteristics of the multibody structure accurately.

Two models with different gripper finger lengths were investigated. As shown in Fig. 3, the large case is with a total length of 1800 mm, and the small case is with a total length of 300 mm. Besides, each segment length of the large one is 600 mm while the small one is 100 mm. The small spring is with 6 mm diameters and 20 rolls, the large one is with 20 mm and 20 rolls. The selection of the length and structure of the large case is inspired by elephant trunk. The structural cross-section exhibits a triangular configuration, akin to the anatomical structure of an elephant's trunk, characterized by notable flexibility and a significant load-bearing capacity. Moreover, this particular geometric configuration enables substantial bending deformations along specific orientations essential for these design specifications [37]. The small case is scaled by the large model to facilitate the conduction of experiments and validations easily. Also, the modeling dimensions of all other components of the small gripper are scaled accordingly.