Skip to content

Computational design

Besides modeling, the field of computational inverse design can also benefit from the use of FEM models. Building up the \class{Fem} class, the objective is to find a topological structure of a continuum system based on a desired deformations or compliance. One widely adopted method is the Solid Isotropic Material with Penalization (SIMP) approach, which is a commonly used material interpolation technique in topology optimization \cite{Bendsoe2003}. In the SIMP method, each finite element \(e \in \{1,2,...,n_e \}\) is assigned a continuous density variable \(\rho_e \in (0,1]\), which serves as an indicator of the material distribution within the mesh. If \(\rho_e = 1\), the element is considered solid, while if \(\rho_e = 0\), the element is considered void. This assignment of density variables enables the modification of the strain energy density in \(\Psi\): \(\renewcommand{\elastic}{_\textrm{e}}\) \(\renewcommand{\grav}{_\textrm{g}}\) \(\renewcommand{\p}{\partial}\) \(\renewcommand{\R}{\mathbb{R}}\) \(\renewcommand{\vec}[1]{\mathbf{#1}}\) \(\renewcommand{\mat}[1]{\mathbf{#1}}\) \(\renewcommand{\FB}{\mathbf{F}}\) \(\renewcommand{\LB}{\mathbf{L}}\) \(\renewcommand{\fB}{\mathbf{f}}\) \(\renewcommand{\Hm}{\mathcal{H}}\) \(\renewcommand{\xB}{\mathbf{x}}\) \(\renewcommand{\AB}{\mathbf{A}}\) \(\renewcommand{\x}{\mathbf{x}}\) \(\renewcommand{\dxB}{\dot{\mathbf{x}}}\) \(\renewcommand{\dt}{\delta t}\) \(\renewcommand{\ddxB}{\ddot{\mathbf{x}}}\)

\[ \tilde{\Psi}_{e} = \left[ \varepsilon + (1-\varepsilon){\rho_e}^p \right] \Psi, \]

where \(0 < \varepsilon \ll 1\) a lower bound on the densities, and \(p > 1\) a penalty factor for penalizing intermediate densities during the optimization process. By collecting the density values \(\vec{\rho} = \textrm{col}\{\rho_1,\rho_2,...,\rho_{N_e}\}\), the inverse design problem can be formulated in terms of two unknowns: the displacement field \(\xB\) and the density field \(\vec{\rho}\). Consequently, the computational design problem for general soft material structures can be expressed as a nonlinear topology optimization problem of the following form:

\[ \begin{align} \underset{\vec{\rho}}{\textrm{minimize}}\quad & {\Phi = -\beta_1\, \LB^\top\!\x \;+\; \beta_2 \fB\elastic^\top \! \fB_{\textrm{u}} } \\ \textrm{subject to} \quad& \tilde{\vec{r}}(\xB,\vec{\rho})=0, \\[0.25em] & \vec{v}^\top \vec{\rho}\le v^\star, \\[0.25em] & \vec{\rho}\in \mathcal{P}, \label{eq:C5:topo_optimization} \end{align} \]

where \(\LB\) a sparse unit-vector composed of nonzero entries for the degrees-of-freedom corresponding to the desired morphology of the soft robot, \(\vec{v}\) the element volumes, \(v*\) the desired volume infill, \(\mathcal{P} = \left\{ \vec{\rho} \in \R^{n_e} ~|~ 0 < \rho_i \le 1\right\}\) admissible set for the design variables, and \(\beta_1\) and \(\beta_2\) are positive scalars that can be adjusted to vary the optimization problem, with \(\beta_1 \ll \beta_2\) resulting in compliance minimization and \(\beta_1 \gg \beta_2\) leading to a compliant mechanism. To solve the optimization problem in \eqref{eq:C5:topo_optimization}, we utilize the Method of Mixed Asymptotes (MMA) proposed by Svanberg12. Earlier work on this computational design approach was presented in Caasenbrood et al. 3, and the chapter is referred to this work for the analytic gradients required for the MMA solver.

The optimization routine in the Sorotoki framework is incorporated into the Fem class and can be invoked by utilizing the command fem.optimize('type'), where 'type' represents the optimization problem at hand. For minimizing compliance, the cost function is self-adjoint \cite{Bendsoe2003}, hence objective function and constraints are linear operators. However, when dealing with compliant mechanisms, it is necessary to specify the selection vector \(\LB\), which can be defined using the fem.addOutput(id) command. The value of id represents the nodal indices of interest, which can be identified using the fem.Mesh.findNode functionality.

Example: optimization of linear elastic beam

Image title

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
msh = Mesh(sRectangle(4, 10),'Quads',[40,100]);
msh = msh.generate();

fem = Fem(msh,'SpatialFilterRadius',0.3);

fem = fem.addMaterial(NeoHookean);
fem = fem.addSupport('nw',[1, 1]);
fem = fem.addSupport('sw',[1, 1]);
fem = fem.addLoad('leftmid',[1,0]);

fem.options.isNonlinear = false;
fem.options.LineStyle = 'none';
fem.options.Display = @plt;

fem = fem.optimize;

function plt(Fem)
    cla;
    showInfillFem(Fem);
end


  1. Krister Svanberg. The method of moving asymptotes\ifmmode —\else —\fi a new method for structural optimization. International Journal for Numerical Methods in Engineering, 24(2):359–373, 1987. doi:10.1002/nme.1620240207

  2. Krister Svanberg. Mma and gcmma-two methods for nonlinear optimization. vol, 1:1–15, 2007. 

  3. Brandon Caasenbrood, Alexander Pogromsky, and Henk Nijmeijer. Dynamic modeling of hyper-elastic soft robots using spatial curves. IFAC Proceedings Volumes (IFAC-PapersOnline), 2020.