1 Introduction↩︎

Soft actuators are compliant, flexible structures that generate motion or force through the deformation of soft materials such as elastomers, hydrogels, or shape-memory polymers. These actuators have demonstrated significant potential across a broad range of applications, including assisting surgeons during minimally invasive procedures ([1]), manipulating delicate objects ([2]) and navigating hazardous or constrained environments ([3], [4]). Various actuation mechanisms have been developed for these systems, including pneumatic and hydraulic pressure ([5], [6], [7]), electric fields ([8], [9]), chemical reactions ([10]) and magnetic fields ([11]). Among these approaches, fluidic actuation has emerged as one of the most widely adopted methods due to its structural simplicity, ability to achieve large deformations and straightforward fabrication processes.

Traditionally, the design of soft actuators is based on heuristics, relying on engineering intuition and parameter space exploration ([6], [12]). Naturally, such methods are potentially incapable of producing designs that achieve specific target behaviors. In this regard, design methods based on gradient-based optimization are much more efficient, whereby the optimized material distribution is provided under given objective and constraint functions. Gradient-based optimization in the form of shape or topology optimization has proven its use for designing pneumatic soft actuators using implementations based on finite strain hyperelasticity ([13], [14] and [15]) and linear elasticity ([16], [17] and [18]). However, the aforementioned works are limited to shape morphing in two dimensions, which naturally limits the methods capabilities of finding truly novel designs with enhanced capabilities. Exception to this rule is the work of [19], who integrates gradient-based topology optimization with a material point method to design locomotive soft robots, and [20] who proposes a shape optimization scheme based on B-splines to design soft actuators with desired deformation behavior.

Shape optimization is an excellent candidate for the design of pneumatic soft actuators, since the pressurized internal cavities of the actuators naturally lend themselves to shape morphing. Indeed, topology optimization should in theory be able to explore a more vast design space, but it is not entirely clear how to simultaneously control each independent air chamber’s actuation pressure. Also, the gas or fluid pressure must somehow be transferred to the solid body, which is nontrivial when the design boundaries are a priori unknown, and entails additional modeling of e.g. incompressible fluid regions ([21], [22]) or porohyperelasticity ([13], [23]). A shape optimization framework has previously been developed to design three dimensional soft actuators by [20]. Their approach is based on tracking the evolving surfaces of the hyperelastic body during optimization using B-spline surfaces, and they regulate the surface quality by imposing geometric penalty constraints. Their computational workflow is not unified, but subdivided such that different pieces of software handle geometry, finite element modeling and optimization. This comes with a high computational cost, wherefore they must rely on solving a sequence of subproblems that are constructed within a trust region in which the cost and constraint functions are approximated by their first-order Taylor polynomials.

The present work proposes a unified computational framework capable of designing shape optimized three dimensional soft actuators, taking into account both material and geometrical nonlinearities. The computational framework is based on a nearly incompressible hyperelastic finite deformation (FE) implementation on unstructured meshes that utilizes the Portable and Extendable Toolkit for Scientific Computing (PETSc) for efficiency ([24]). The shape optimization is simple in the sense that it is parameter-free; it directly utilizes the nodal coordinates in the finite element mesh as design variables. In this approach to shape optimization, originally proposed by [25], the mesh quality is regulated via a non-linear PDE filter ([26], [27]). The design is updated using the gradient-based method: Method of Moving Asymototes (MMA, [28]), and the gradients of the cost and constraint functions are obtained from an adjoint sensitivity analysis. We demonstrate the effectiveness of the proposed framework by designing actuators capable of performing nontrivial deformations upon actuation, such as object grasping, contraction under pressurization and staggered multimodal deformation modes. Finally, we validate the computational framework by fabricating the optimized designs using a casting approach and experimentally characterize their performance.

2 Hyperelastic formulation↩︎

The response of the hyperelastic body is captured using a Lagrangian kinematic description in which \(\Omega\) is the undeformed configuration with material points \(\text{\boldsymbol{X}}\in\Omega\). A pressure load deforms \(\Omega\) into the current configuration \(\Omega_c=\text{\boldsymbol{\varphi}}(\Omega)\), with spatial points \(\text{\boldsymbol{x}} = \text{\boldsymbol{\varphi}}(\text{\boldsymbol{X}}) = \text{\boldsymbol{X}} + \text{\boldsymbol{u}}(\text{\boldsymbol{X}})\), where \(\text{\boldsymbol{\varphi}}\) is the smooth deformation and \(\text{\boldsymbol{u}}\) is the corresponding displacement. The local deformation is described by the deformation gradient \(\text{\boldsymbol{F}} = \text{\boldsymbol{\nabla}} \text{\boldsymbol{\varphi}} = \textsf{\boldsymbol{1}} + \text{\boldsymbol{\nabla}}\text{\boldsymbol{u}}\), Jacobian \(J(\text{\boldsymbol{F}}) = \text{det}(\text{\boldsymbol{F}})\), right Cauchy-Green deformation tensor \(\text{\boldsymbol{C}}(\text{\boldsymbol{F}}) = \text{\boldsymbol{F}}^T\text{\boldsymbol{F}}\) and the Green-Lagrangian strain \(\text{\boldsymbol{E}}(\text{\boldsymbol{C}}) =\frac{1}{2}\left(\text{\boldsymbol{C}} - \textsf{\boldsymbol{1}}\right)\), where \(\text{\boldsymbol{1}}\) is the second order identity tensor.

2.1 Material model↩︎

We cast the soft actuator from an elastomer that is assumed to be isotropic, hyperelastic and nearly incompressible, wherefore it is modeled using the neo-Hookean strain energy density \[{W}(\text{\boldsymbol{F}}) = \Psi(J(\text{\boldsymbol{F}}), \text{\boldsymbol{C}}(\text{\boldsymbol{F}}))=\frac{1}{2}G\left(J^{-2/3}\text{tr}(\text{\boldsymbol{C}}) -3\right)+ \frac{1}{2}K(J- 1)^2, \label{c1}\tag{1}\] where, in the limit of infinitesimal strain, \(K = \frac{E}{3(1-2v)}\) and \(G = \frac{E}{2(1+\nu)}\) denote the bulk and shear modulii, respectively, and \(E\) is the Young’s modulus. To obtain a nearly incompressible behavior, we let the Poisson’s ratio \(\nu = 0.49\), such that \(J(\text{\boldsymbol{X}})\approx 1\), \(\forall \text{\boldsymbol{X}}\in\Omega\). To resolve the issue of nonphysical oscillations in pressure when the stress is computed directly from the strain energy function for nearly incompressible materials, we use the mixed formulation proposed by [29]. In this formulation, we introduce an independent pressure field \(\tilde{p}\) such that we can define a modified strain energy function \[\widetilde{W}(\text{\boldsymbol{F}}, \tilde{p}) = W(\text{\boldsymbol{F}}) - \frac{1}{2K} \left(p-\tilde{p}\right)^2, \label{Wmod}\tag{2}\] where the hydrostatic pressure, \(p\), is derived from the original constitutive relation, i.e. Eq. 1 as \[p = -\frac{\partial \Psi}{\partial J}=- K(J -1). \label{p}\tag{3}\] For later purposes, we note that the second Piola-Kirchhoff stress tensor is \(\widetilde{\text{\boldsymbol{S}}} = \text{\boldsymbol{F}}^{-1} \frac{\partial \widetilde{W}}{\partial \text{\boldsymbol{F}}} = 2\frac{\partial \widetilde{W}}{\partial \text{\boldsymbol{C}}}\).

2.2 Equilibrium equation↩︎

The equilibrium configuration of the elastomer is defined by those sufficiently smooth fields \(\text{\boldsymbol{u}}\) and \(\tilde{p}\) which, for all admissible virtual fields \(\delta \text{\boldsymbol{u}}\) and \(\delta \tilde{p}\), satisfy \[\displaystyle\delta\Pi(\text{\boldsymbol{u}},\tilde{p};\delta\text{\boldsymbol{u}},\delta\tilde{p}) = \displaystyle\int_\Omega \widetilde{\text{\boldsymbol{S}}}: \delta \text{\boldsymbol{E}} \, dV \displaystyle+ \int_\Omega \frac{1}{K}\left(p - \tilde{p} \right)\delta \tilde{p} \, dV+\int_{\partial\Omega_{c}^{p}}\hat{p} \text{\boldsymbol{n}}_c\cdot\delta\text{\boldsymbol{u}} \, dS =0, \label{deltaW01}\tag{4}\] where \(\delta\text{\boldsymbol{E}} =\frac{1}{2}\left((\text{\boldsymbol{\nabla}}\delta\text{\boldsymbol{u}})^T\text{\boldsymbol{F}} + \text{\boldsymbol{F}}^T \text{\boldsymbol{\nabla}}\delta\text{\boldsymbol{u}}\right)\) is the virtual Lagrangian strain. In the above, the deformed boundary \(\partial\Omega_c\) with the unit normal \(\text{\boldsymbol{n}}_c\), consists of two complementary surfaces, \(\partial\Omega_{c}^{p}\) and \(\partial\Omega^{u}=\partial\Omega_{c}^{u}\), over which the pressure load \(\hat{p}\text{\boldsymbol{n}}_c\) and the null displacement \(\text{\boldsymbol{u}} = \text{\boldsymbol{0}}\), are prescribed, respectively. It is emphasized that \(\hat{p}\) is the relative pressure compared to the atmospheric pressure.

3 Shape optimization↩︎

Our shape optimization is defined similar to the aforementioned hyperelastic problem in the sense that it morphs the initial design \(\Omega_o\) into the shape optimized design \(\Omega\) via a smooth vector field \(\text{\boldsymbol{\psi}}\), which takes points \(\text{\boldsymbol{X}}_o\in\Omega_o\) to points \(\text{\boldsymbol{X}}\in\Omega\) via \(\text{\boldsymbol{X}} = \text{\boldsymbol{\varphi}}_\psi(\text{\boldsymbol{X}}_o) = \text{\boldsymbol{X}}_o + \text{\boldsymbol{\psi}}(\text{\boldsymbol{X}}_o)\). The field \(\text{\boldsymbol{\psi}}\) is driven by the shape optimization design variable vector field \(\text{\boldsymbol{d}}\) via a PDE filter defined by the potential ([25], [30], [27]) \[\Pi_\psi(\text{\boldsymbol{\psi}};\text{\boldsymbol{d}}) = \displaystyle\frac{1}{2}\int_{\Omega_o} W_\psi(\text{\boldsymbol{F}}_\psi) \, dV + \frac{1}{2\vert\Omega_o\vert}\int_{\Omega_o} \vert\vert\text{\boldsymbol{\psi}}-\text{\boldsymbol{d}}\vert\vert^2 \, dV, \label{PDEpot01}\tag{5}\] where \(W_\psi\) is a fictitious strain energy density function. We also introduce the shape field deformation gradient \(\text{\boldsymbol{F}}_\psi = \text{\boldsymbol{\nabla}} \text{\boldsymbol{\varphi}}_\psi= \textsf{\boldsymbol{1}} + \text{\boldsymbol{\nabla}}\text{\boldsymbol{\psi}}\). The smooth shape field \(\text{\boldsymbol{\psi}}\) is found by invoking stationarity in Eq. 5 , i.e. \[\delta\Pi_\psi(\text{\boldsymbol{\psi}};\text{\boldsymbol{d}},\delta\text{\boldsymbol{\psi}}) = \displaystyle\int_{\Omega_o} \text{\boldsymbol{S}}_\psi : \delta\text{\boldsymbol{E}}_{\psi} \, dV + \frac{1}{\vert\Omega_o\vert}\int_{\Omega_o} \left(\text{\boldsymbol{\psi}}-\text{\boldsymbol{d}}\right) \cdot \delta\text{\boldsymbol{\psi}} \, dV = 0, \label{PDEpot02}\tag{6}\] which should be fulfilled for all admissible virtual fields \(\delta\text{\boldsymbol{\psi}}\). We prescribe the null displacement \(\text{\boldsymbol{\psi}} = \text{\boldsymbol{0}}\) over \(\partial\Omega_o^\psi\).

The framework for mapping the initial domain \(\Omega_o\) to the shape optimized domain \(\Omega\) to the deformed domain \(\Omega_c\) is illustrated in Fig. 1.

3.1 Fictitious material model↩︎

The fictitious strain energy density function \(W_\psi\) should be chosen such that it adequately smooths the shape field \(\text{\boldsymbol{\psi}}\). In the literature, several possible energy density functions have been proposed, see also [30]. In this work, we follow [27] and use the neo-Hookean like energy density \[\begin{array}{ll} \displaystyle W_\psi(\text{\boldsymbol{F}}_\psi) &=\displaystyle\Psi_\psi(J_\psi(\text{\boldsymbol{F}}_\psi), {\text{\boldsymbol{C}}}_{\psi}(\text{\boldsymbol{F}}_\psi)) \\[10pt] &= \displaystyle\frac{1}{2}K_\psi\left(\frac{1}{k_\psi}\left(e^{k_\psi\left(J_\psi-1\right)}-1\right) - \text{ln}(J_\psi) \right) + \frac{1}{2}G_\psi\left(J_\psi^{-2/3}\text{tr}({\text{\boldsymbol{C}}}_\psi) - 3 \right), \end{array} \label{Wshape}\tag{7}\] where \(\text{\boldsymbol{C}}_\psi(\text{\boldsymbol{F}}_\psi) = \text{\boldsymbol{F}}^T_\psi\text{\boldsymbol{F}}_\psi\), \(J_\psi(\text{\boldsymbol{F}}_\psi) = \text{det}(\text{\boldsymbol{F}}_\psi)\) and \(K_\psi\), \(G_\psi\) and \(k_\psi=10\) are filter parameters. The above fictitious energy density was proposed by [27] and is tailored to avoid self penetration and excessive volume change and promote designs with good mesh quality and smooth design boundaries.

4 FE-formulations↩︎

We discretize the equilibrium equations in Eqs. 4 and 6 using finite elements. The discretization procedure of the hyperelastic problem is described in Section 4.1, whereas the shape optimization problem is considered in Section 4.2.

4.1 Hyperelastic PDE↩︎

We discretize Eq. 4 using 10-1 tetrahedral mixed \(\text{\boldsymbol{u}}\)-\(p\) finite elements. In our total Lagrangian FE approach, we interpolate the displacement field in each finite element \(e\) using the shape functions (see e.g. [31]), i.e. \(\text{\boldsymbol{u}}(\text{\boldsymbol{X}}) \approx \textsf{\boldsymbol{N}}(\text{\boldsymbol{X}})\textsf{\boldsymbol{u}}^e\), where \(\textsf{\boldsymbol{u}}^e\) contains the nodal displacement components, whereas the independent pressure field is interpolated as element piece-wise uniform, i.e. \(\tilde{p}(\text{\boldsymbol{X}}) \approx \tilde{\textsf{p}}^e\) for \(\text{\boldsymbol{X}}\in\Omega^e\) where \(\tilde{\textsf{p}}^e\) is the element pressure. The virtual strains are discretized as \(\delta\textsf{\boldsymbol{E}} = \textsf{\boldsymbol{B}} \delta\textsf{\boldsymbol{u}}^e\), where the explicit format for the discrete strain operator \(\textsf{\boldsymbol{B}}\) appear in e.g. [31].

Introducing the aforementioned discretization and using the arbitrariness of \(\delta\textsf{\boldsymbol{u}}\) and \(\delta\tilde{\textsf{\boldsymbol{p}}}\) in Eq. 4 yields the nonlinear residual equation \[\textsf{\boldsymbol{r}}_{{a}}(\textsf{\boldsymbol{a}}, \lambda) = \begin{bmatrix} \displaystyle\textsf{\boldsymbol{r}}_{{u}} \\[10pt] \displaystyle\textsf{\boldsymbol{r}}_p \end{bmatrix} = \begin{bmatrix} \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\Bigg(\displaystyle\int_{\Omega^e} \textsf{\boldsymbol{B}}^T\widetilde{\textsf{\boldsymbol{S}}} \, dV +\lambda\int_{\partial\Omega^{p,e}_{c}} \hat{p} \textsf{\boldsymbol{N}}^T \text{\boldsymbol{n}}_c \, dS \Bigg) \\[10pt] \displaystyle\mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega^e}\frac{1}{K}(p-\tilde{p}) \, dV \end{bmatrix} = \textsf{\boldsymbol{0}}, \label{res}\tag{8}\] where \(\mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\) is the finite element assembly operator and we introduce \(\textsf{\boldsymbol{a}}=[\textsf{\boldsymbol{u}},\;\tilde{\textsf{\boldsymbol{p}}}]^T\) and the load scaling parameter \(\lambda\in[0,1]\). During loading, the equilibrium response in Eq. 8 can exhibit non-linear phenomena such as snapping and bifurcations. Therefore we use an arclength method, wherein we explicitly enforce the auxiliary hypersurface constraint \[r_c = \text{\boldsymbol{\Delta}}\textsf{\boldsymbol{a}}^T\text{\boldsymbol{\Delta}}\textsf{\boldsymbol{a}} - \alpha^2 = 0, \label{arclength}\tag{9}\] in every loadstep ([32]). In the above, \(\alpha > 0\) is a parameter that controls the magnitude of the load increment and \(\text{\boldsymbol{\Delta}}\textsf{\boldsymbol{a}} = \textsf{\boldsymbol{a}}^{(n)} - \textsf{\boldsymbol{a}}^{(n-1)}\), where the superscript \(n\) represents the load step. In every load step, we linearize Eq. 8 about the current iterate \(\textsf{\boldsymbol{a}}\) and \(\lambda\) \[\textsf{\boldsymbol{r}}_{{a}}(\textsf{\boldsymbol{a}} + d\textsf{\boldsymbol{a}}, \lambda + d\lambda)\approx \textsf{\boldsymbol{r}}_{{a}}(\textsf{\boldsymbol{a}}, \lambda) + \frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial \textsf{\boldsymbol{a}}}d\textsf{\boldsymbol{a}} + \frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial \lambda}d\lambda = \textsf{\boldsymbol{0}}, \label{Newt1}\tag{10}\] and solve for the displacement and load increments \(d\textsf{\boldsymbol{a}}\) and \(d\lambda\) via \[\frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial \textsf{\boldsymbol{a}}}d\textsf{\boldsymbol{a}} + \frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial \lambda}d\lambda = \textsf{\boldsymbol{K}} d\textsf{\boldsymbol{a}} + \textsf{\boldsymbol{P}}d\lambda = -\textsf{\boldsymbol{r}}_{{a}}, \quad \text{where} \quad \textsf{\boldsymbol{K}} = \begin{bmatrix} \textsf{\boldsymbol{K}}_{{uu}} & \textsf{\boldsymbol{K}}_{{u}p} \\[10pt] \textsf{\boldsymbol{K}}_{p{u}} & \textsf{\boldsymbol{K}}_{pp} \end{bmatrix} \quad \text{and} \quad \textsf{\boldsymbol{P}} = \begin{bmatrix} \textsf{\boldsymbol{P}}_{u} \\[10pt] \textsf{\boldsymbol{0}} \end{bmatrix}. \label{Newt2}\tag{11}\] Ultimately, we solve two linear systems \(\textsf{\boldsymbol{K}} d\textsf{\boldsymbol{a}}_r = - \textsf{\boldsymbol{r}}_a\) and \(\textsf{\boldsymbol{K}} d\textsf{\boldsymbol{a}}_f = -\textsf{\boldsymbol{P}}\) and let \(d\textsf{\boldsymbol{a}} = d\textsf{\boldsymbol{a}}_r + d\lambda d\textsf{\boldsymbol{a}}_f\). We then obtain \(d\lambda\) via algebraic manipulations of Eq. 9 ; for details we refer to [32]. In Eq. 11 we introduce \[\begin{array}{ll} \displaystyle\textsf{\boldsymbol{K}}_{{uu}} = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega^e} \left(\textsf{\boldsymbol{B}}^T\widetilde{\textsf{\boldsymbol{D}}}\textsf{\boldsymbol{B}} + \textsf{\boldsymbol{G}}^T\widetilde{\textsf{\boldsymbol{Y}}}\textsf{\boldsymbol{G}} \right) \, dV + \lambda \textsf{\boldsymbol{K}}^p, \\[10pt] \displaystyle\textsf{\boldsymbol{K}}_{{u}p} = \displaystyle\textsf{\boldsymbol{K}}_{p{u}}^{T} = -\mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega^e} \textsf{\boldsymbol{B}}^TJ\textsf{\boldsymbol{C}}^{-1}\, dV, \\[10pt] \displaystyle\textsf{\boldsymbol{K}}_{pp} = -\mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega^e} \frac{1}{K}\, dV,\\[10pt] \displaystyle\textsf{\boldsymbol{P}}_{u} = \int_{\partial\Omega^{p,e}_{c}} \hat{p} \textsf{\boldsymbol{N}}^T \text{\boldsymbol{n}}_c \, dS, \end{array} \label{LinMat}\tag{12}\] where \(\widetilde{\textsf{\boldsymbol{D}}}\) is the material tangent tensor \(\mathbb{D} = 4 \frac{\partial^2 \widetilde{W}}{\partial \text{\boldsymbol{C}}\partial \text{\boldsymbol{C}}}\) expressed in Voigt notation and the explicit formats of \(\widetilde{\textsf{\boldsymbol{Y}}}\) and \(\textsf{\boldsymbol{G}}\) appear in e.g. [31]. The pressure load tangent contribution \(\textsf{\boldsymbol{K}}^p\) is defined in the Appendix. We emphasize that our tangent stiffness matrix is symmetric, since the pressure is uniform and is applied to a sufficiently constrained boundary surface ([33] and [34]).

4.1.1 Static condensation↩︎

Ultimately, we do not solve the full systems in Eq. 11 explicitly, but instead utilize the piece-wise discontinuous interpolation of the pressure field such that \(d\tilde{{\textsf{p}}}^e\) is statically condensed out on the element level. To illustrate this procedure, consider the elementwise expansion of \(\textsf{\boldsymbol{K}} d\textsf{\boldsymbol{a}}_r = - \textsf{\boldsymbol{r}}_{a}\), i.e. \[\begin{bmatrix} \textsf{\boldsymbol{K}}_{{uu}}^e & \textsf{\boldsymbol{K}}_{{u}p}^e \\[10pt] \textsf{\boldsymbol{K}}_{p{u}}^e & {K}_{pp}^e \end{bmatrix}\begin{bmatrix} d\textsf{\boldsymbol{u}}^e\\[10pt] d\tilde{\textsf{p}}^e \end{bmatrix} = -\begin{bmatrix} \textsf{\boldsymbol{r}}_{{u}}^e\\[10pt] {r}_p^e \end{bmatrix}. \label{statCond0}\tag{13}\] The second row of Eq. 13 reveals \[d\tilde{{\textsf{p}}}^e = -\left(\textsf{K}_{pp}^e\right)^{-1}\left(\textsf{r}_p^e + \textsf{\boldsymbol{K}}_{p{u}}^e d\textsf{\boldsymbol{u}}^e\right), \label{statCond1}\tag{14}\] which, when inserted in the first row of Eq. 13 , yields \[\textsf{\boldsymbol{K}}^e d\textsf{\boldsymbol{u}} = \left(\textsf{\boldsymbol{K}}_{{uu}}^e - \textsf{\boldsymbol{K}}_{{u}p}^e \left(\textsf{K}_{pp}^e\right)^{-1}\textsf{\boldsymbol{K}}_{p{u}}^e\right)d\textsf{\boldsymbol{u}}^e = -\textsf{\boldsymbol{r}}_{{u}}^e + \textsf{\boldsymbol{K}}_{{u}p}^e\left(\textsf{K}_{pp}^e\right)^{-1}{r}_p^e=\textsf{\boldsymbol{r}}^e. \label{statCond2}\tag{15}\] where \(\textsf{\boldsymbol{K}}^e\) and \(\textsf{\boldsymbol{r}}^e\) are the condensed element stiffness matrix and residual vector, respectively. We solve \(\textsf{\boldsymbol{K}} d\textsf{\boldsymbol{a}}_p = -\textsf{\boldsymbol{P}}\) using the same methodology.

4.2 Shape filter PDE↩︎

We use 10 node tetrahedral single field finite elements to solve Eq. 6 . Again, the vector fields are interpolated via the shape functions \(\textsf{\boldsymbol{N}}\), such that \(\text{\boldsymbol{d}}(\text{\boldsymbol{X}}_o) \approx \textsf{\boldsymbol{N}}(\text{\boldsymbol{X}}_o)\textsf{\boldsymbol{d}}^e\) and \(\require{upgreek} \text{\boldsymbol{\psi}} (\text{\boldsymbol{X}}_o)\approx \textsf{\boldsymbol{N}}(\text{\boldsymbol{X}}_o)\text{\boldsymbol{\uppsi}}^e\), where \(\textsf{\boldsymbol{d}}^e\) and \(\require{upgreek} \text{\boldsymbol{\uppsi}}^e\) are the element nodal shape displacement and filtered shape displacement vectors. Using the arbitrariness of \(\delta\text{\boldsymbol{\psi}}\), the discretized version of Eq. 6 requires \(\require{upgreek} \text{\boldsymbol{\uppsi}}\) to satisfy \[\require{upgreek} \textsf{\boldsymbol{r}}_\psi(\text{\boldsymbol{\uppsi}};\textsf{\boldsymbol{d}}) = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\left(\displaystyle\int_{\Omega_o^e} \textsf{\boldsymbol{B}}^T\textsf{\boldsymbol{S}}_\psi \, dV - \displaystyle\int_{\partial\Omega_o^{d,e}} \textsf{\boldsymbol{N}}^T\left(\text{\boldsymbol{d}}-\text{\boldsymbol{\psi}}\right) \, dS \right)= \textsf{\boldsymbol{0}}, \label{PDE04}\tag{16}\] To solve Eq. 16 , we utilize a traditional Newton’s method. For more details we refer to [27].

5 Optimization problem↩︎

We use our shape optimization framework to design soft actuators with bespoke deformation patterns. To this end, our multiobjective shape optimization problem reads \[(\mathbb{SO}) \;\begin{cases} \underset{\textsf{\boldsymbol{d}}}{\text{min}} \;\underset{i \in[1,n_g]}{\text{max}} \;\;\{g_i\}, \\[10pt] \text{s.t} \quad \begin{cases} f \leq 0, \quad \\ \textsf{d}_k\in[\underline{\textsf{d}}, \overline{\textsf{d}}], \quad k\in[1, 3n_n], \end{cases} \end{cases} \label{opt}\tag{17}\] where \(n_g\) is the number of objectives, \(n_n\) is the number of finite element nodes and the design shape variables are constrained to be in the range of \(\underline{\textsf{d}} \leq \textsf{d}_j \leq \overline{\textsf{d}}\). The range of the box constraints is a delicate choice made by the user. In the above, \(g_i\) are objective functions and \(f\) is an inequality constraint. In our numerical experiments, we follow [27] and use rather tight bounds, but increase the design freedom by performing several spaced updates of the mesh coordinates such that \(\text{\boldsymbol{X}}_o \leftarrow \text{\boldsymbol{X}}_o + \text{\boldsymbol{\psi}}\) and \(\textsf{\boldsymbol{d}} \leftarrow \text{\boldsymbol{0}}\). Our experience is that this approach produces designs with better mesh quality compared to those designs obtained using large bounds without mesh coordinate updates. We emphasize that the equilibrium equality constraints \(\textsf{\boldsymbol{r}}_{{a}} = \textsf{\boldsymbol{0}}\) and \(\textsf{\boldsymbol{r}}_{\psi} = \textsf{\boldsymbol{0}}\) are explicitly enforced, i.e. we solve the optimization problem in a staggered fashion. The gradient-based nonlinear programming method MMA is used to solve the optimization problem defined in Eq. 17 , and we utilize its min-max formulation ([28]).

We seek designs of soft actuators which exhibits desired deformation behaviors over a surface \(\partial\Omega^d\subset \partial\Omega\) at target pressures. To this end, each objective is defined as \[g_i = \frac{w_i}{\vert\partial\Omega^d\vert}\int_{\partial\Omega^d} \text{\boldsymbol{u}}^{i} \cdot \text{\boldsymbol{v}}^i \, dS, \label{g0}\tag{18}\] where \(w_i\) are weights, \(\vert\partial\Omega^d\vert = \int_{\partial\Omega^d} \, dS\) is the area, \(\text{\boldsymbol{u}}^{i}\) is the displacement field and \(\text{\boldsymbol{v}}^i\) is the vector direction field at target point \(i\). In Eq. 17 , we also introduce the inequality constraints \[f = \frac{1}{\vert\partial\Omega^d\vert}\int_{\partial\Omega^d} \vert\text{\boldsymbol{u}}^{j} \cdot \text{\boldsymbol{s}}^j \vert^2 \, dS - \epsilon_u \leq 0, \label{constraint}\tag{19}\] which constrains the deformation at target point \(j\) in the direction \(\text{\boldsymbol{s}}^j\) to be smaller than the relaxation tolerance \(\epsilon_u \geq 0\). By solving the optimization problem as posed in Eq. 17 , we promote designs which exhibit tailored deformation patterns.

5.1 Sensitivity analysis↩︎

We compute the gradients of a function \(\require{upgreek} \tilde{g}(\textsf{\boldsymbol{d}}) =g\left({\text{\boldsymbol{\uppsi}}}(\textsf{\boldsymbol{d}}),\textsf{\boldsymbol{a}}({\text{\boldsymbol{\uppsi}}}(\textsf{\boldsymbol{d}}))\right)\) using the adjoint method. In this method, we augment \(\tilde{g}\) with the equality constraints (\(\textsf{\boldsymbol{r}}_{{a}} = \textsf{\boldsymbol{0}}\) and \(\textsf{\boldsymbol{r}}_{{{\psi}}} = \textsf{\boldsymbol{0}}\)) via the associated Lagrange multipliers \(\text{\boldsymbol{\mu}}_{{a}}\) and \(\text{\boldsymbol{\mu}}_{{{\psi}}}\) to obtain the identity \[\bar{g} := \tilde{g} - \text{\boldsymbol{\mu}}_{{a}}^T\textsf{\boldsymbol{r}}_{{a}} - \text{\boldsymbol{\mu}}_{{{\psi}}}^T\textsf{\boldsymbol{r}}_{{{\psi}}}. \label{sens01}\tag{20}\] Next, we differentiate Eq. 20 with respect to \(\textsf{\boldsymbol{d}}\) and rearrange, to obtain \[\require{upgreek} \displaystyle\frac{d\bar{g}}{d \textsf{\boldsymbol{d}}} = - \text{\boldsymbol{\mu}}_{{{\psi}}}^{T}\frac{\partial \textsf{\boldsymbol{r}}_{{{\psi}}}}{\partial \textsf{\boldsymbol{d}}}+ \Bigg[\frac{\partial g}{\partial \text{\boldsymbol{\uppsi}}} - \text{\boldsymbol{\mu}}_{{a}}^T\frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial {\text{\boldsymbol{\uppsi}}}} - \text{\boldsymbol{\mu}}_{{{\psi}}}^{T}\frac{\partial \textsf{\boldsymbol{r}}_{{{\psi}}}}{\partial {\text{\boldsymbol{\uppsi}}}} +\left(\frac{\partial g}{\partial \textsf{\boldsymbol{a}}}- \text{\boldsymbol{\mu}}_{{a}}^T\frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial \textsf{\boldsymbol{a}}} \right)\frac{\partial \textsf{\boldsymbol{a}}}{\partial {{\text{\boldsymbol{\uppsi}}}}} \Bigg]\frac{\partial {\text{\boldsymbol{\uppsi}}}}{\partial \textsf{\boldsymbol{d}}}. \label{sens03}\tag{21}\] To annihilate the implicit derivatives \(\require{upgreek} \frac{\partial \textsf{\boldsymbol{a}}}{\partial {{\text{\boldsymbol{\uppsi}}}}}\) in Eq. 21 , we solve the linear system \[\left(\frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial \textsf{\boldsymbol{a}}}\right)^T \text{\boldsymbol{\mu}}_{{a}} = \left(\frac{\partial g}{\partial \textsf{\boldsymbol{a}}}\right)^T, \label{adj2}\tag{22}\] to obtain \(\text{\boldsymbol{\mu}}_{{a}}\). We emphasize that we again use static condensation when solving the above linear system. Inserting Eq. 22 in Eq. 21 gives \[\require{upgreek} \displaystyle\frac{d\bar{g}}{d \textsf{\boldsymbol{d}}} = \displaystyle- \text{\boldsymbol{\mu}}_{{{\psi}}}^{T}\frac{\partial \textsf{\boldsymbol{r}}_{{{\psi}}}}{\partial \textsf{\boldsymbol{d}}} + \Bigg[\frac{\partial g}{\partial {{\text{\boldsymbol{\uppsi}}}}} - \text{\boldsymbol{\mu}}_{{a}}^T\frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial {{\text{\boldsymbol{\uppsi}}}}} - \text{\boldsymbol{\mu}}_{{\psi}}^{T}\frac{\partial \textsf{\boldsymbol{r}}_{{{\psi}}}}{\partial {{\text{\boldsymbol{\uppsi}}}}} \Bigg]\frac{\partial {{\text{\boldsymbol{\uppsi}}}}}{\partial \textsf{\boldsymbol{d}}}. \label{sens04}\tag{23}\] The implicit derivative \(\require{upgreek} \frac{\partial {{\text{\boldsymbol{\uppsi}}}}}{\partial \textsf{\boldsymbol{d}}}\) is annihilated by solving the adjoint problem \[\require{upgreek} \left(\frac{\partial \textsf{\boldsymbol{r}}_{{{\psi}}}}{\partial {\text{\boldsymbol{\uppsi}}}}\right)^T \text{\boldsymbol{\mu}}_{{{\psi}}} = \left(\frac{\partial g}{\partial {\text{\boldsymbol{\uppsi}}}}\right)^T - \left(\frac{\partial \textsf{\boldsymbol{r}}_{{a}}}{\partial {\text{\boldsymbol{\uppsi}}}}\right)^T\text{\boldsymbol{\mu}}_{{a}}, \label{adj3}\tag{24}\] for \(\text{\boldsymbol{\mu}}_{{{\psi}}}\). Finally, the sensitivity expression reduces to \[\displaystyle\frac{d \bar{g}}{d \textsf{\boldsymbol{d}}} = - \text{\boldsymbol{\mu}}_{{{\psi}}}^{T}\frac{\partial \textsf{\boldsymbol{r}}_{{{\psi}}}}{\partial \textsf{\boldsymbol{d}}}. \label{sens06}\tag{25}\] We emphasize that since we utilize an arclength method to solve the equilibrium equations, care must be taken to ensure that we always compute the sensitivity at set target pressures. To this end, we monitor the pressure during loading, and when we surpass a target pressure, we return to the previous load step and perform a standard force-controlled Newton step to exactly end up at the set target pressure. More details for the sensitivity computations are provided in the Appendix.

6 Experiments↩︎

To validate the proposed shape optimization framework, we fabricate and test several of the inversely designed soft actuators. The fabrication and experimental testing procedure are described in detail in the following two sections.

6.1 Fabrication↩︎

We cast our soft actuators using polyvinyl siloxane (PVS, Zhermack Elite Double 32). For each tested design, we fabricate and test two actuators to ensure repeatability. The molds are design in Blender based on the surface representations (STL) of our optimized designs, which we extract using ParaView. We adopt a compound molding technique in which the molds consists of rigid positive and flexible negative parts. The rigid positive molds are 3D printed using the Formlabs Form 3B with 0.05 mm layer thickness using the Model V3 resin. After the print is finished, the molds are washed in isopropanol in the Formlabs Wash V1 for 10 minutes to remove excess resin. Subsequently, any support structure is manually removed, after which the part is washed for another 5 minutes. The molds are then cured at room-temperature for at least 24 hours.

The soft actuator are casted via a two-step molding process, as illustrated in Fig. 2. First, inspired by the work of [35], we cast the flexible negative mold from the PVS. 3D printed molds are used to cast the complementary soft negative mold of the actuator’s closed cavity geometry using PVS. Before casting, we insert two metal rods in the bottom of the mold which will act as alignment pins when mounting, and we apply a layer of mold release agent (Mann Ease Release 200) on the positive mold (Fig. 2A). Then we mix the two-part, base and catalyst, PVS elastomer (THINKY mixer ARE-310) at 2000 RPM for 30 seconds, and subsequently at 2200 RPM for 30 seconds, to ensure uniformity in its composition and to degas. We then pour the mixed PVS inside the mold (Fig. 2B), which is clamped and placed at room temperature for 20 minutes before opening it. After demolding (Fig. 2C), the negative flexible mold (with alignment pins, see Fig. 2D) is mounted on the rigid positive base plate of the actual actuator mold (Fig. 2E), thereby defining the external geometry of the actuator to form the final casting assembly (Fig. 2E). We do not bake the flexible negative mold. We acknowledge that the compliance of the PVS mold introduces small manufacturing errors, but have found that for intricate cavity shapes, the de-molding might otherwise cause severe rupture of the actuator, wherefore a soft negative mold proves beneficial.

The second step is to cast the actual actuators. First, we apply a layer of mold release agent (Mann Ease Release 200) to both the positive and negative molds. Then we mix PVS elastomer using the same technique as explained above. The mixed PVS is then poured inside the mold (Fig. 2F), which is clamped and placed at room temperature for 20 minutes. The demolding process consist of first removing the positive mold (Fig. 2G), and then carefully peeling the actuator from the negative mold (Fig. 2H). The actuators (Fig. 2I) are subsequently baked at \(60\) \(°\)C for 24 hours in an oven to stabilize the mechanical properties.

The actuators are finally mounted on a laser-cut (Universal Laser Systems PLS6.150D with a 75W C02 lase) acrylic support plate, which is engraved to increase roughness and promote bonding. To do this, we first apply a primer (Aron Alpha PP Primer) to the base of the actuator and to the acrylic support plate to promote adhesion. Then we use Locktite SuperGlue to glue the actuator to the acrylic plate. The glue is left to cure for 24 hours.

6.2 Testing↩︎

The actuators are tested using a syringe pump (Harvard Apparatus 33DS), pressure sensor (15 PSI Ashcroft GV Pressure Transducer), data acquisition device (DAQ, Saleae Logic pro 8), tubing, and mount for the actuator (Fig. 3 (a)). Using the syringe pump, we perform volume-controlled loading of air with a flow rate of \(10\) ml\(/\)min. Both the pressure signal and volume input are recorded as a function of time using the DAQ. During the tests, actuator deformation was recorded using a digital camera (Nikon Z 6II). Black dot markers were placed on the actuator surface, and their positions were tracked in the recorded videos using the Kanade-Lucas-Tomasi (KLT) feature-tracking algorithm implemented in MATLAB’s Computer Vision Toolbox.

To compare the numerically predicted and experimentally measured pressure–volume curves, it is necessary to account for both the compliance of the experimental setup and the compressibility of air. The compliance of the setup, which originates from the tubing and syringe, is quantified by measuring the pressure–volume response of the system without the actuator. Specifically, 10 mL of air is infused at a constant flow rate of 10 mL/min, and the corresponding pressure response is recorded. The measured pressure–volume curve, along with a quadratic fit, is presented in Fig. 3 (b). The quadratic fit can be extrapolated to provide the volume correction due to the system compliance, \(\Delta V_{compliance}\), such that \[\Delta V_{corr} = \Delta V_{syringe} - \Delta V_{compliance}, \label{compliance}\tag{26}\] where \(\Delta V_{syringe}\) is the amount of volume dispensed by the syringe and \(\Delta V_{corr}\) is the system compliance corrected volume change.

The compressibility of air is accounted for under two assumptions: (i) the flow rate is sufficiently low for the process to be considered isothermal, and (ii) the enclosed air obeys the ideal gas law. Under these assumptions, the volume change of the pressurized cavity of the actuator, \(\Delta V^p\), can be expressed as (see, e.g., [7], [36] for details) \[\Delta V^p = \Delta V_{corr} - \frac{\hat{p}}{\hat{p}+\hat{p}_0} V^p, \label{dV}\tag{27}\] where \(\hat{p}_0\) denotes the atmospheric pressure and \(V^p\) is the initial cavity volume. Throughout this study, the volume change of the pressurized cavity in our experiments is computed using Eq. 27 .

7 Results↩︎

We demonstrate our shape optimization framework by morphing the initial geometry depicted in Fig. 4 to achieve a range of target deformations. The initial geometry consists of a closed lower (blue) cylindrical cavity with an open (gray) cylindrical cavity positioned on top. The cylinder is clamped at the lower surface \(\partial\Omega^u\). Pressurized gas is injected into the lower closed cylindrical cavity defined by the surface \(\partial\Omega^p\). We target various deformation modes of the open surface \(\partial\Omega^d\) of the top open cylinder by employing the objectives and constraints defined in Eqs. 18 and 19 . The target deformation modes include extension, contraction, and grasping through closure of the top cylindrical cavity. The geometrical parameters defining the initial geometry are summarized in Table I.

To solve general 3D design problem, we require parallel computing in the form of message passing (MPI) routines and linear algebra tools from PETSc ([24]). We generate the unstructured finite element grids using Gmsh ([37]), and the discretizations are distributed using ParMETIS ([38]) and are managed by DMPlex structures in PETSc. The linear systems arising from the discretized PDE:s are solved using the direct solver MUMPS ([39]), which provides the Cholesky factorizations. All numerical solutions are obtained using the COSMOS cluster at Lund University and are run on one node consisting of two AMD EPYC 7413 processors (48 CPU cores @ 3.6 GHz) and 256 GB of ram.

In all numerical examples, we set the MMA move limit to \(0.025\). The optimization is terminated after 350 design iterations, since we after this observe minimum changes in the design and objective. The optimizer is restricted from morphing the shape of the actuator at the attachment surface, \(\partial \Omega^u\), and at the top cylindrical cavity surface, \(\partial \Omega^d\), i.e. \(\partial\Omega_o^\psi = \partial \Omega^u \cup \partial \Omega^d\). The finite element mesh consists of \(109604\) quadratic tetrahedral elements, resulting in \(n_n = 187611\) nodes. The material properties of the PVS are \(E = 1.2\) MPa and \(\nu = 0.49\) ([36]). We omit the constraint \(f\) in the optimization problem definition in Eq. 17 , unless stated otherwise.

We emphasize that the choices of \(K_\psi\), \(G_\psi\), \(\overline{\textsf{d}}\) and \(\underline{\textsf{d}}\) are problem dependent and may require some numerical investigation. For example, large bounds on the box constraints on the design variables \(\overline{\textsf{d}}\) and \(\underline{\textsf{d}}\) and small \(K_\psi\) and \(G_\psi\) modulii enables large shape modifications, with the risk of producing designs with tangled meshes and overlapping surfaces. Decreasing the bounds and increasing the modulii mitigates this effect, at the expense of limiting the extent of the shape modifications. Indeed, a trade-off between allowing larger shape changes and ensuring a reasonable mesh quality exists. In the following examples, we fix \(K_\psi = 1 \times 10^{-5}\) mm\(^{-1}\) and \(\overline{\textsf{d}} = -\underline{\textsf{d}} = 1.2\) mm, but alter \(G_\psi\) in Eq. 7 to accommodate the problem specific mesh quality.

7.1 Inverse design of a gripping actuator↩︎

In the first example, we design an actuator capable of grasping an object by maximizing the inward radial displacement of the cylindrical surface \(\partial\Omega^d\) shown in Fig. 4 upon application of a pressure \(\hat{p} = 15\) kPa to the closed cavity. To this end, we minimize the objective function \(g_1\) defined in Eq. 18 with \(\text{\boldsymbol{v}}^1 = (0, -\cos(\theta), -\sin(\theta))\) and \(w_1 = -1\), where \(\theta\) denotes the angular coordinate in the \(e_y\)-\(e_z\) plane (Fig. 6). For this problem, we consider a single objective, which is computed at the maximum pressure, i.e. at a load scaling factor of \(\lambda = 1\) (see Eq. 8 ). We update the mesh coordinates every 50 design iterations, for a total of three updates, and use a filter parameter of \(G_\psi = 1\times10^{-2}\) mm\(^{-1}\) in Eq. 7 .

The history of the objective function evolution is plotted in Fig. 9 (a), where also snapshots of the intermediate designs are included. We observe that the objective function monotonically decreases, with drastic drops appearing in conjunction with every mesh coordinate update, as expected. The optimized design and its cast realization is depicted in Fig. 7 (b) in its undeformed configuration at \(\hat{p}=0\) kPa. The design is depicted in its deformed configuration at \(\hat{p}=15\) kPa in Fig. 7 (c). These snapshots indicate that, in the optimal design, two arch-shaped bulges beneath the top cylinder, \(\partial\Omega^d\), expand outward during pressurization and induce the observed elliptical deformation pattern of \(\partial\Omega^d\). The top views of the open cavity shown in Fig. 7 (d) clearly indicate that the initially circular opening transforms into an elongated elliptical shape, thereby facilitating . Fig. 7 (b)-(c) also includes sliced views of the actuator in the \(e_z\)-\(e_x\) plane, further illustrating the deformation mechanism.

To validate the computational framework, we experimentally characterize the response of the identified optimized geometry. We perform motion tracking of the two black markers highlighted in Fig. 7 (d) and compare the experimental and simulated displacements over time (Fig. 7 (e)), see also the movie provided as the supplementary material. Although the overall behavior agrees well (see the deformation patterns in the deformed configuration in Fig. 7 (d)), we note slight discrepancies between the numerical and experimental motion tracking results. We believe that these errors stem from manufacturing errors, such as mold misalignment and air bubbles, that invariably appear in the manufacturing process. Also, camera misalignment can play a role in the discrepancies, especially since the magnitude of the displacements is relatively small. To further characterize the gripping actuator, we compare the numerically computed and experimentally measured pressure-volume relation. The result is depicted in Fig. 7 (f), where we plot the applied pressure \(\hat{p}\) versus volume change of the pressurized cavity¹, i.e. \(\Delta V^p = V^p_c - V^p\). We observe a nearly linear pressure-volume relation, and a good correlation between the numerical and experimental results (the experiments were conducted three times each for the two fabricated actuators). Overall, the numerical and experimental results show good agreement, indicating that the proposed numerical framework accurately captures the behavior of the designed actuators.

Finally, Fig. 8 demonstrates the capabilities of the gripping actuator through a series of grasping experiments. The actuator successfully grasps and lifts several objects, including a thread spool (4.7 g), a whiteboard marker (12.5 g), and a metallic ball (27.6 g). All objects have diameters smaller than 20 mm, requiring the actuator to actively grasp rather than simply support them. We note that the actuator is pressurized manually using a syringe and that the lifting motion is performed by hand. Despite these simple experimental conditions, the actuator is capable of reliably grasping and lifting objects of varying shapes and weights.

7.2 Inverse design of linear actuators↩︎

In the next example, we design actuators that maximizes the output displacement parallel to the \(e_x\)-direction upon pressurization to \(\hat{p}=15\) kPa, thereby acting as linear actuators. We perform two optimizations to design 1) an extending actuator by setting \(\text{\boldsymbol{v}}^1=(1,0,0)\) in Eq. 18 and 2) a contracting actuator by setting \(\text{\boldsymbol{v}}^1=(-1,0,0)\). For both cases, we employ a single objective function and \(w_1 = -1\), evaluated at a load scaling factor of \(\lambda=1\). Furthermore, the mesh coordinates are updated every 50 design iterations, resulting in a total of three mesh updates. The filter parameter is set to \(G_\psi=3\times10^{-2}\) mm\(^{-1}\) in Eq. 7 .

Results for the extending actuator are shown in Fig.  9. Again, we observe a monotonically decreasing objective function \(g_1\) (Fig. 9 (a)). The identified optimal design and its sliced view in the \(e_z\)-\(e_x\) plane are depicted in Fig. 9 (b) in the undeformed configuration at \(\hat{p} = 0\) kPa and in Fig. 9 (c) in the deformed configuration at \(\hat{p} = 15\) kPa. Unsurprisingly, the optimizer morphs the shape of the actuator to form bulges, which results in expansion in the positive \(e_x\)-direction upon pressurization. In Fig. 9 (d), we plot the objective function \(g_1\), which corresponds to the area average surface integral of the displacement over the surface \(\partial\Omega^d\), versus the applied pressure \(\hat{p}\). We find that \(g_1\) monotonically decreases with increasing pressure and that the largest expansion magnitude occurs at the maximum pressure. We emphasize that the objective \(g_1\) is negative, as we have a minimization optimization problem.

Results for the contracting actuator are shown in Fig. 10. The evolution of the objective function is presented in Fig. 10 (a), together with snapshots of the intermediate designs. The final optimized design is shown in Fig. 10 (b). The optimization morphs the geometry such that a cross-section in the \(e_y\)–\(e_z\) plane exhibits four curved walls arranged in a cross-shaped pattern. Upon pressurization, these walls expand outward within the plane, causing the upper cylindrical section to contract in the negative \(e_x\) direction, as illustrated in Fig. 10 (c).

Interestingly, the maximum magnitude of the desired displacement occurs at a pressure below the target value of \(\hat{p}=15\) kPa. This behavior is evident in Fig. 10 (d), which shows the objective function \(g_1\) as a function of the applied pressure \(\hat{p}\). We attribute this phenomenon to a ballooning effect: above a certain pressure level, the initially curved walls shown in Fig. 10 (b) progressively approach the nearly square deformation pattern depicted in Fig. 10 (c). As a result, further pressurization leads to displacement in the opposite of the intended direction. The optimizer could potentially improve the design by increasing the curvature of these walls, thereby enhancing the contraction mechanism. However, the design-variable bounds and filtering constraints limit this possibility. Attempts to relax these constraints resulted in designs exhibiting self-intersection. Although the operating pressure could be reduced to avoid the ballooning effect, we choose to retain the example in its current form to illustrate this behavior.

7.3 Inverse design of actuators that support sequential deformation modes↩︎

Finally, we consider actuators capable of sequential deformation under an applied pressure of \(\hat{p}=25\) kPa. Specifically, we design two actuators: one that first grasps an object and subsequently extends in the \(e_x\) direction, and another that first grasps an object and subsequently contracts in the \(e_x\) direction. For both optimization problems, we employ two objective functions (i.e., \(n_g=2\) in Eq. 17 ). For the first objective, we set \(\text{\boldsymbol{v}}^1=(0,-\cos(\theta),-\sin(\theta))\) in Eq. 18 and evaluate it at a load scaling factor of \(\lambda^1=0.25\). For the second objective, we set \(\text{\boldsymbol{v}}^2=(1,0,0)\) for axial expansion and \(\text{\boldsymbol{v}}^2=(-1,0,0)\) for axial contraction, and evaluate it at a load scaling factor of \(\lambda^2=1.0\). To promote the desired sequential behavior, we set \(\text{\boldsymbol{s}}^1=(1,0,0)\) and \(\varepsilon_u=1\) in the constraint defined by Eq. 19 . This constraint mitigates deformation in the \(e_x\) direction at the first target point, where the actuator is optimized only for the grasping motion. Although the min-max formulation is used within the MMA optimizer, the two objectives must still be weighted differently to achieve the desired behavior. In this example, we set \(w_1=-18\) and \(w_2=-1\). The filter parameter is set to \(G_\psi=3\times10^{-2}\) mm\(^{-1}\) in Eq. 7 , and the mesh coordinates are updated every 50 design iterations, resulting in a total of four mesh updates.

Results for the actuator that that first grasps and subsequently elongates are shown in Fig. 11. The evolution of the objective functions \(g_1\) and \(g_2\), as well as the constraint \(f\), along with snapshots of the intermediate designs is plotted versus the design iterations in Fig. 11 (a). The identified optimal design and a sliced view of it is depicted in Fig. 11 (b) in the undeformed configuration at \(\hat{p} = 0\) kPa. The deformed design at the first target pressure \(\hat{p} = 6.25\) kPa is depicted in Fig. 11 (c), and at the second target pressure \(\hat{p} = 25\) kPa in Fig. 11 (d). It is seen that at the first target pressure, the design barely displaces in the \(e_x\) directions, as expected by the introduction of the constraint; it focuses on performing the gripping deformation as seen in the top-view in Fig. 11 (c). At the second target pressure, we instead observe the desired displacement in the positive \(e_x\) direction, see Fig. 11 (d). We curiously note that a larger gripping deformation occurs at the maximum pressure \(\hat{p} = 25\) kPa, compared to that of the first, lower, target pressure \(\hat{p} = 6.25\) kPa where this deformation mode is actually maximized. Based on this observation, we assume that the actuator design that maximizes the gripping deformation at the maximum pressure, also maximizes this deformation mode at the lower pressure. In Fig. 11 (e), we plot of the objectives \(g_1\) and \(g_2\) versus the applied pressure \(\hat{p}\). Herein, we confirm the above observations; that the magnitude of \(g_2\) is small around \(\hat{p} = 6.25\) kPa where we introduce the constraint, and that the magnitude of \(g_1\) is larger around \(\hat{p} = 25\) kPa than around \(\hat{p} = 6.25\) kPa.

Finally, in Fig. 12 we report numerical and experimental results for an actuator that first grasps and subsequently contracts. The evolution of the objective functions \(g_1\) and \(g_1\), as well as the constraint \(f\), is plotted versus the design iterations in Fig. 12 (a), where also snapshots of the intermediate designs are included. A side-view and a top-view of the optimized design and its cast realization are depicted in Fig. 12 (b) and Fig. 12 (c) in their undeformed configuration at \(\hat{p}=0\) kPa, and deformed configurations at target pressures \(\hat{p}=3.75\) kPa and \(\hat{p}=15\) kPa. Again, we perform motion tracking of the points highlighted in Fig. 12 (b)-(c) and compare the experimental and simulated displacements over time. The results are depicted in Fig. 12 (d)-(e) and in the movie provided in the supplementary material, where we plot the displacements trajectories during pressurization. A glance at Fig. 12 (b)-(c) confirms that the overall experimental deformation patterns agrees well with the simulation. However, we again note slight discrepancies between the numerical and experimental displacements trajectories in Fig. 12 (d)-(e). To further characterize the actuator, we compare the numerically computed and experimentally measured pressure-volume relation in Fig. 12 (f), which show good agreement. Similarly to the contracting actuator, we observe a non-linear pressure-volume relation.

8 Conclusions↩︎

In this work, we design soft actuators using a three-dimensional shape optimization framework. The simulation-based framework is founded on hyperelasticity, nonlinear kinematics, and mixed finite elements, and is implemented in PETSc. The numerical model is shown to agree well with experimental measurements. While minor local deviations are observed in the deformation response, the global behavior exhibits excellent agreement. The optimized designs demonstrate tailored staggered deformations in both simulations and experiments.

Future work will focus on incorporating additional manufacturing constraints into the optimization framework. Demolding has proven to be a critical step in the fabrication process, and curvature constraints may improve manufacturability [40]. Furthermore, the robustness of the shape optimization framework could be enhanced by mitigating surface self-intersections through contact-aware optimization techniques [41]. Another important direction is the incorporation of contact mechanics into the optimization framework to enable the direct optimization of grasping performance. In this context, the third-medium contact method [42], [43] represents a promising approach for modeling contact interactions.

9 Acknowledgments↩︎

This work was performed under the auspices of the Swedish research council (grant nbr. 2024-00172). The numerical computations were enabled by resources provided by LUNARC, The Centre for Scientific and Technical Computing at Lund University. Last but not least, the authors would like to thank Alex Zhang and Leon Kamp for help with the experimental setup.

Appendix↩︎

Evaluation of pressure load↩︎

To evaluate the external pressure load, we introduce the isoparametric mapping \(\text{\boldsymbol{x}} = \textsf{\boldsymbol{N}}_\xi(\xi_1, \xi_2)\textsf{\boldsymbol{x}}^e\) for all \(\text{\boldsymbol{x}}\in\partial\Omega_c\), where \(\textsf{\boldsymbol{N}}_\xi\) are the element shape functions corresponding to a surface element, \(\textsf{\boldsymbol{x}}^e\) are the deformed element coordinates and \(\xi_1\) and \(\xi_2\) are the isoparametric coordinates. We introduce the deformed current domain \(\Omega_c\) with coordinates \(\text{\boldsymbol{x}}\in\Omega_c\) and boundary \(\partial\Omega_c\) with normal \(\text{\boldsymbol{n}}_c\), defined as \[\text{\boldsymbol{n}}_c = \frac{\frac{\partial \text{\boldsymbol{x}}}{\partial \xi_1}\times\frac{\partial \text{\boldsymbol{x}}}{\partial \xi_2}}{\left\lVert\frac{\partial \text{\boldsymbol{x}}}{\partial \xi_1}\times\frac{\partial \text{\boldsymbol{x}}}{\partial \xi_2}\right\rVert_2} = \frac{\hat{\text{\boldsymbol{n}}}_c}{\left\lVert\hat{\text{\boldsymbol{n}}}_c\right\rVert_2}. \label{normal}\tag{28}\] The element external load due to the pressure load is therefore \[\textsf{\boldsymbol{F}}_{p} = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\partial\Omega_c} p \textsf{\boldsymbol{N}}_\xi^T \text{\boldsymbol{n}}_c \, dS = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega_\xi} p \textsf{\boldsymbol{N}}_\xi^T \frac{\hat{\text{\boldsymbol{n}}}_c}{\left\lVert\hat{\text{\boldsymbol{n}}}_c\right\rVert_2} \left\lVert\hat{\text{\boldsymbol{n}}}_c\right\rVert_2 \, d\xi_1 d\xi_2 = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega_\xi} p \textsf{\boldsymbol{N}}_\xi^T \hat{\text{\boldsymbol{n}}}_c \, d\xi_1 d\xi_2, \label{isoPressureLoad}\tag{29}\] where \(\Omega_\xi\) is the isoparametric domain, and we used that \(dS = \left\lVert\hat{\text{\boldsymbol{n}}}\right\rVert_2 \, d\xi_1 d\xi_2\). In the above, \(p = 0.2\) is the total applied pressure. The stiffness matrix contribution from the pressure load is \[\textsf{\boldsymbol{K}}_p = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\frac{\partial {\textsf{\boldsymbol{F}}_p}}{\partial \textsf{\boldsymbol{u}}^e} = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega_\xi} p \textsf{\boldsymbol{N}}_\xi^T \frac{\partial \hat{\text{\boldsymbol{n}}}_c}{\partial \textsf{\boldsymbol{u}}^e} \, d\xi_1 d\xi_2, \label{dPdshape}\tag{30}\] where the components of \(\frac{\partial \hat{\text{\boldsymbol{n}}}_c}{\partial \textsf{\boldsymbol{u}}^e}\) are \[\begin{array}{ll} \displaystyle\frac{\partial \hat{n}_i}{\partial \textsf{u}_p^\beta} &\displaystyle= \varepsilon_{ijk} \left(\frac{\partial }{\partial \textsf{u}_p^\beta}\left(\frac{\partial x_j}{\partial \xi_1}\right)\frac{\partial x_k}{\partial \xi_2} + \frac{\partial x_j}{\partial \xi_1}\frac{\partial }{\partial \textsf{u}_p^\beta}\left(\frac{\partial x_k}{\partial \xi_2}\right) \right) \\[15pt] &\displaystyle= \varepsilon_{ijk} \left(\frac{\partial N^\alpha}{\partial \xi_1} \frac{\partial \textsf{u}_j^\alpha}{\partial \textsf{u}_p^\beta}\frac{\partial x_k}{\partial \xi_2} + \frac{\partial x_j}{\partial \xi_1}\frac{\partial N^\alpha}{\partial \xi_2} \frac{\partial \textsf{u}_k^\alpha}{\partial \textsf{u}_p^\beta} \right)\\[15pt] &\displaystyle= \varepsilon_{ipj} \frac{\partial N^\beta}{\partial \xi_1} \frac{\partial x_j}{\partial \xi_2} + \varepsilon_{ijp} \frac{\partial x_j}{\partial \xi_1}\frac{\partial N^\beta}{\partial \xi_2}\\[15pt] &\displaystyle= \varepsilon_{ipj} \left(\frac{\partial N^\beta}{\partial \xi_1} \frac{\partial x_j}{\partial \xi_2} - \frac{\partial x_j}{\partial \xi_1}\frac{\partial N^\beta}{\partial \xi_2}\right). \end{array} \label{dnhatdshape}\tag{31}\] and subscript and superscript indices denote dimensional and nodal indices, respectively.

Shape sensitivity of pressure load↩︎

The pressure force sensitivity is \[\frac{\partial \text{\boldsymbol{\mu}}_a^T\textsf{\boldsymbol{F}}_{p}}{\partial \text{\boldsymbol{\psi}}^e} = \mathop{ \ooalign{ \hfil\displaystyle\sum\cr \hfil\textstyle\sum\cr}}\int_{\Omega_\xi} p \left(\text{\boldsymbol{\mu}}_a^e\right)^T\textsf{\boldsymbol{N}}_\xi^T \frac{\partial \hat{\text{\boldsymbol{n}}}_c}{\partial \text{\boldsymbol{\psi}}^e} \, d\xi_1 d\xi_2, \label{shapePressureLoad}\tag{32}\] where \(\frac{\partial \hat{\text{\boldsymbol{n}}}_c}{\partial \text{\boldsymbol{\psi}}^e}\) coincides with 31 . The shape sensitivity of the internal force appears in [14].

Shape sensitivity of objective function↩︎

The shape sensitivity of the objective function is \[\frac{\partial g_i}{\partial \text{\boldsymbol{\psi}}^e} = -\frac{1}{\vert\partial\Omega^d\vert^2} \frac{\partial \vert\partial\Omega^d\vert}{\partial \text{\boldsymbol{\psi}}^e} \int_{\partial\Omega^d} \text{\boldsymbol{u}} \cdot \text{\boldsymbol{v}} \, dS+ \frac{1}{\vert\partial\Omega^d\vert} \int_{\partial\Omega^d} \text{\boldsymbol{u}} \cdot \frac{\partial \text{\boldsymbol{v}}}{\partial \text{\boldsymbol{\psi}}^e} \, dS + \frac{1}{\vert\partial\Omega^d\vert} \int_{\partial\Omega^d} \text{\boldsymbol{u}} \cdot \text{\boldsymbol{v}} \, \frac{\partial dS}{\partial \text{\boldsymbol{\psi}}^e}, \label{g0sens}\tag{33}\] where \(\frac{\partial \text{\boldsymbol{v}}}{\partial \text{\boldsymbol{\psi}}^e}\) is problem specific and \[\frac{\partial \vert\partial\Omega^d\vert}{\partial \text{\boldsymbol{\psi}}^e} = \int_{\partial\Omega^d}\, \frac{\partial dS}{\partial \text{\boldsymbol{\psi}}^e}. \label{dSsens}\tag{34}\] In the above, we evaluate \[\frac{\partial dS}{\partial \text{\boldsymbol{\psi}}^e} = \frac{\partial \left\lVert\hat{\text{\boldsymbol{n}}}\right\rVert_2}{\partial \text{\boldsymbol{\psi}}^e} \, d\xi_1 d\xi_2 = \frac{1}{\left\lVert\hat{\text{\boldsymbol{n}}}\right\rVert_2} \frac{\partial \hat{\text{\boldsymbol{n}}}}{\partial \text{\boldsymbol{\psi}}^e} \cdot \hat{\text{\boldsymbol{n}}} \, d\xi_1 d\xi_2 = \frac{\partial \hat{\text{\boldsymbol{n}}}}{\partial \text{\boldsymbol{\psi}}^e} \cdot \text{\boldsymbol{n}} \, d\xi_1 d\xi_2, \label{dSsens2}\tag{35}\] where the components of \(\frac{\partial \hat{\text{\boldsymbol{n}}}}{\partial \text{\boldsymbol{\psi}}^e}\) are \[\begin{array}{ll} \displaystyle\frac{\partial \hat{n}_i}{\partial \psi_p^\beta} &\displaystyle= \varepsilon_{ijk} \left(\frac{\partial }{\partial \psi_p^\beta}\left(\frac{\partial X_j}{\partial \xi_1}\right)\frac{\partial X_k}{\partial \xi_2} + \frac{\partial X_j}{\partial \xi_1}\frac{\partial }{\partial \psi_p^\beta}\left(\frac{\partial X_k}{\partial \xi_2}\right) \right) \\[15pt] &\displaystyle= \varepsilon_{ijk} \left(\frac{\partial N^\alpha}{\partial \xi_1} \frac{\partial \psi_j^\alpha}{\partial \psi_p^\beta}\frac{\partial X_k}{\partial \xi_2} + \frac{\partial X_j}{\partial \xi_1}\frac{\partial N^\alpha}{\partial \xi_2} \frac{\partial \psi_k^\alpha}{\partial \psi_p^\beta} \right)\\[15pt] &\displaystyle= \varepsilon_{ipj} \frac{\partial N^\beta}{\partial \xi_1} \frac{\partial X_j}{\partial \xi_2} + \varepsilon_{ijp} \frac{\partial X_j}{\partial \xi_1}\frac{\partial N^\beta}{\partial \xi_2}\\[15pt] &\displaystyle= \varepsilon_{ipj} \left(\frac{\partial N^\beta}{\partial \xi_1} \frac{\partial X_j}{\partial \xi_2} - \frac{\partial X_j}{\partial \xi_1}\frac{\partial N^\beta}{\partial \xi_2}\right). \end{array} \label{dnhatdshape2}\tag{36}\]

References↩︎

---

1. We compute the current cavity volume \(V^p_c\) using the divergence theorem, i.e. \(V^p_c = \int_{\Omega^p} dV = \frac{1}{3}\int_{\Omega^p} \text{\boldsymbol{\nabla}}_c\cdot \text{\boldsymbol{x}} \, dV = \frac{1}{3}\int_{\partial\Omega^p} \text{\boldsymbol{x}} \cdot \text{\boldsymbol{n}}_c \, dS\).↩︎