Soft materials, particularly rubber-like elastomers, have attracted significant attention due to their ability to undergo large deformations while maintaining structural integrity. Their exceptional mechanical properties-high flexibility, adaptability, and resilience-make them well-suited for applications requiring safe and efficient interaction with dynamic environments. These materials have been widely applied in various fields, including soft robotics [1,2], tissue modeling [3], and flexible electronics [4], enabling innovations in areas that demand deformability and multifunctionality. Among these applications, soft robotics [2,5] has emerged as a transformative field, enabling the design of robotic systems that surpass the limitations of traditional rigid structures. Recent advancements in soft robotics have demonstrated significant potential across various domains, including biomedical applications [6], wearable assistive devices [7], and marine exploration [8]. By harnessing the intrinsic deformability of soft materials, soft robots offer enhanced adaptability and improved safety in human-centered and unstructured environments.
To accurately capture their nonlinear deformation behavior, various hyperelastic material models [9], such as the Neo-Hookean, Mooney-Rivlin, and Ogden models [10], have been developed. Unlike linear elastic materials, hyperelastic materials are governed by constitutive relations derived from a strain energy density function. These models use different numbers of material parameters to describe specific mechanical responses [11], enabling precise representation of diverse hyperelastic behaviors. For instance, the Neo-Hookean model uses a single parameter, making it suitable for simple rubber-like materials, while the Mooney-Rivlin and Ogden models incorporate multiple parameters, allowing them to represent a wider range of hyperelastic responses.
While traditional soft robot designs have primarily relied on intuition and bio-inspired principles [6,12], the field's advancement has increasingly demanded a systematic approach to achieve optimal designs. Topology optimization (TO), originally developed for the aerospace industry [13], offers such a methodology and has been widely adopted across various physical domains, including thermal [14,15], magnetic [16,17], and electrical structures [18].
Topology optimization employs a variety of methods, broadly categorized into density-based approaches like SIMP [19,20] and its variations, evolutionary strategies [21] such as BESO [22], phase field methods [23,24], component-based parametric methods like Moving Morphable Components (MMC) [25]/Moving Morphable Voids(MMV) [26], and boundary-based techniques, most notably level set methods [19,20]. While SIMP is popular for its simplicity, its element-wise density representa-tion can lead to numerical artifacts like checkerboarding and indistinct boundaries. Phase Field Methods offers a smooth transition between material phases. MMC/MMV relies on geometrical components, which is good to include manufacturing constraints. Level set methods, however, offer a distinct advantage by implicitly defining structural boundaries, resulting in smoother, more readily manufacturable designs. Consequently, this research focuses on the level set method, specifically to address its currently limited application in the domain of hyperelastic topology optimization.
Level set methods define material boundaries implicitly, enabling smooth and precise boundary evolution. Despite its advantages, its application to hyperelastic topology optimization remains limited, leaving a gap in computational design for soft robotics and compliant mechanisms.
Several researchers have made progress in hyperelastic topology optimization. Choi and Duan [27] and Kim et al. [28] developed continuum-based shape sensitivity analysis for hyperelastic materials, utilizing the material derivative concept and shape design parameterization to compute shape sensitivities. However, parameterized shape optimization imposes inherent limitations on topological changes, restricting the emergence of new structural features.
To address this limitation, Cho applied level set methods to nonlinear topology optimization [29,30]. However, this approach neglected high-order terms in the material derivatives, particularly those related to higher-order displacement gradients, which may introduce inaccuracies in the calculated velocity field, leading to imprecise boundary evolution.
Later, Feifei Chen [31] improved level set-based topology optimization for hyperelastic materials by incorporating geometric and material nonlinearities. Retaining higher-order displacement gradient terms enhanced accuracy and convergence, while using linear topology optimization as an initial step improved numerical stability. However, Chen's approach did not directly optimize hyperelastic structures, and primarily focused on mean compliance optimization. Given the broader applicability of hyperelastic materials in compliant mechanisms and soft robotics, a more generalized topology optimization framework is still needed.
Building on existing shape sensitivity results, this paper rederives the shape sensitivity analysis using the adjoint method and incorporates it into the level set framework. While extensive research has focused on minimizing the mean compliance of cantilever beam in linear designs using the level set method [32,33], its application to hyperelastic materials remains less explored. To assess the effectiveness of the proposed approach, the mean compliance problem is first solved for these benchmark structures. Subsequently, a displacement inverter compliant mechanism [34] is designed using the derived sensitivity analysis, fully leveraging the unique advantages of hyperelastic materials.
The paper is structured as follows: Section 2 presents the modeling of hyperelastic materials, followed by Section 3, which introduces the level set method used in topology optimization. Next, the shape sensitivity analysis is derived. Section 4 presents numerical examples including applying the derived sensitivity analysis to a cantilever beam problem aiming to minimize the mean compliance of the structure (Section 4.1). Finally, Section 4.2 presents the results for the displacement inverter compliant mechanism, demonstrating the full potential of hyperelastic materials in compliant mechanism design.
The mechanical behavior of hyperelastic materials plays a crucial role in topology optimization for soft robotics and compliant mechanisms. Unlike linear elastic materials, hyperelastic materials exhibit large deformations while maintaining elasticity, requiring specialized constitutive models to accurately describe their stress-strain relationships. Several hyperelastic models have been developed, including Mooney-Rivlin, Ogden, and neo-Hookean formulations, each with varying levels of complexity and material parameters.
Among these, the neo-Hookean model is widely used due to its mathematical simplicity and ability to capture essential nonlinear elasticity. In this study, the neo-Hookean material behavior is implemented and simulated using Abaqus/Standard to model large-deformation responses of soft structures. Originally proposed by Treloar (1948), this model describes the strain energy density function in terms of the first invariant of the right Cauchy-Green deformation tensor and the Jacobian determinant of deformation:
where the K bulk modulus and G is the shear modulus from the material. The stress-strain relation is then expressed as:
The tangent stiffness can be expressed as follows: Before proceeding with the derivation of shape sensitivity, it is important to first clarify the fundamental concepts and notation related to hyperelastic materials. This section provides essential background information on hyperelasticity, including key mathematical expressions and symbols used to describe large deformation behavior. By establishing these preliminaries, the following discussions can build upon a clear and consistent foundation.
As shown in Figure 1, the hyperelastic structure deforms into a new shape under surface traction f on boundary Γ. For continuum deformation, any vector on the undeformed structure X is mapped into x by function χ:
The reference configuration (Ω 0 ) and current configuration (Ω) are related by the deformation gradient tensor as introduced:
The variation of deformation gradient tensor F is shown as follow:
Then the material derivative of the deformation gradient tensor Ḟ is:
Where the material derivative of the displacement field is expressed as:
Similarly, the material derivative of the variation of deformation gradient tensor Ḟ is shown as follows:
In this derivation, Green-Lagrange strain tensor E and 2nd Piola Kirchhoff Stress are introduced in the following terms:
Where C is the right Cauchy-Green Tensor. The variation of Green-Lagrange strain tensor E is:
The material derivative of variation of Green-Lagrange strain tensor Ė becomes:
Because S is a symmetric matrix, all sym symbol in Ė and E can be eliminated when they do double dot product with S:
3 Topology Optimization In this section, the level set method is introduced first. Subsequently, we present the formulation of topology optimization for hyperelastic materials using the level set approach. Finally, we discuss the corresponding sensitivity analysis using the adjoint method. The level set method provides an implicit representation of structural geometry using a higher-dimensional function, known as the level-set function Φ. This function defines the structural domain through its iso-surfaces, with the zero level-set contour delineating the material boundary. By intersecting the isosurface with the zero plane, the design domain is naturally partitioned into three regions-material, boundary, and void-based on the sign of Φ, as illustrated in Figure 2. Mathematically, this partition is described by:
In topology optimization, the aim is to evolve the material boundary to reach an optimal configuration. The evolution of this boundary is governed by a velocity field derived from the sensitivity analysis of the objective function in conjunction with the underlying governing equations. Iterative updates to Φ are performed by solving the Hamilton-Jacobi equation, which is a partial differential equation describing the propagation of the level set function over time:
where V n (x) = ẋ is the normal velocity field of the boundary.
This subsection outlines the mathematical formulation for the nonlinear optimization problem of compliant mechanisms. The goal is to balance the objectives of achieving desired deformation and maintaining structural integrity by minimizing compliance. The optimization problem is mathematically stated as:
The variables and parameters involved are defined as follows:
T: Transformation matrix indicating target deformation directions at specific points. S: Second Piola-Kirchhoff stress tensor, representing stresses defined in the reference configuration. ω 1 , ω 2 : Weighting factors satisfying ω 1 + ω 2 = 1. V (Ω): Total volume of the optimized structure. |D|: Prescribed volume fraction.
The nonlinear topology optimization of compliant mechanism inverters can be formulated as a PDE-constrained optimization problem. To solve this constrained optimization problem efficiently, the method of Lagrange multipliers is employed, resulting in the following Lagrangian formulation:
Here, J denotes the original objective function, a(u, v) represents the governing PDE constraint in variational form, and the terms involving λ 1 and λ are the corresponding Lagrange multipliers enforcing the PDE and volume constraints, respectively. Then the material derivative of Lagrangian equation becomes [19,20,35]:
The material derivatives of domain and boundary integrals follow these general expressions: For domain integrals,
For boundary integrals,
Here, the dot notation ( ḟ ) represents material derivatives, while prime notation ( f ′ ) indicates partial derivatives with respect to spatial coordinates. By applying the equation ( 24) and ( 25), the derivatives of integrals encountered in sensitivity analysis can be expressed using either the material derivative of the integrand or its partial spatial derivative. In practice, the derivatives of integral terms can be expressed either as material derivatives (denoted by a dot, e.g., Ḟ, Ė) or spatial partial derivatives (denoted by prime notation, e.g., F ′ , E ′ ). The material derivative representation is typically more straightforward and directly linked to the velocity field V, while the spatial derivative representation often results in more complex expressions. Thus, this study consistently adopts the material derivative notation for clarity and computational convenience.
Then the time material derivative of the objective function:
The body force is not applied, then all terms include b equals zero. The traction force on the right bottom corner is consistent. As a result, all terms include ḟ are eliminated. Because traction force is taken as nodal force, then term κu • fVn equals zero. Finally, the time derivative of the objective function simplifies to:
The time material derivative of the structural energy form:
Substitute the expression of Ė, E and Ė with the expressions above:
Similar to the objective function, some terms are eliminated because of the body force, a derivative of the traction force, and the traction force taken as the nodal force. The time material derivative of the load form:
The time material derivative of volume ration constraint:
Collecting all terms containing u to generate adjoint equation, whose test variable is u and adjoint variable is w. The body force is not considered in this case so the adjoint equation becomes:
By converting the double inner dot into matrix product for simplification, the adjoint equation becomes:
By canceling out the test variable ue and replace the integration by Gauss quadrature, the adjoint equation becomes format Kw e =ω 1 T eω 2 f e to solve for w e value:
After eliminating the terms including u (adjoint equation), and ẇ, considering a(u, ẇ) = l( ẇ), the material derivative of Lagrange equation L becomes:
Similar to reformatting of adjoint equation, the double inner product of Lagrange equation is converted to matrix product and finally the Lagrange equation L is formulated as: By substituting all the results into (35), it yields,
The decent normal velocity field can be derived from the reformulated Lagrange equation simply:
4 Numerical Examples This section presents two numerical examples to demonstrate the effectiveness of the proposed topology optimization framework for hyperelastic materials. The shape sensitivity analysis derived in Section 3 is applied to optimize structural performance. The topology optimization process is implemented in MATLAB, where the design update and sensitivity analysis are carried out iteratively. The Neo-Hookean material properties are set as shear modulus G = 400 kPa and bulk modulus K = 10 MPa.
To validate the proposed method, a benchmark problem-mean compliance minimization of a cantilever beam is investigated. The case assesses the ability of the approach to generate optimal designs considering a mean compliance objective function. Finally, the displacement inverter compliant mechanism is designed to showcase the applicability of hyperelastic topology optimization in compliant mechanism design.
The first numerical example evaluates the accuracy and effectiveness of the proposed shape sensitivity analysis in optimizing the mean compliance of a cantilever beam composed of a hyperelastic material. The objective is to determine the optimal material distribution that minimizes structural compliance while ensuring material efficiency.
The cantilever beam, as illustrated in Figure 3, has dimensions of (2m × 1m) and is subjected to a downward point load of p = (0, -10)N at right bottom corner. The left edge of the structure is fixed. The left edge of the beam is fully fixed, and the thickness is set to 0.1m , assuming a plane stress condition. To discretize the design domain, a 100 × 50 element mesh is employed, and the target volume fraction is set to 50%. The displacement weighting factor is assigned as ω 1 = 0, reducing the objective function to:
The corresponding adjoint equation is formulated as:
The volume constraint remains unchanged throughout the optimization process.
shown in Compared to linear topology optimization, the hyperelastic topology optimization for the mean compliance problem requires more iterations to converge. However, its advantage lies in its ability to directly generate optimized designs from an arbitrary initial topology, ensuring that nonlinear material behavior is captured throughout the optimization process.
The convergence history, depicting the evolution of the objective function and volume ratio, is presented in Figure 5. The objective function exhibits initial fluctuations before stabilizing after approximately 100 iterations, ultimately converging to a lower value, indicating improved structural performance. Meanwhile, the volume ratio gradually decreases and stabilizes near the target 50%, ensuring compliance with the volume constraint. These trends demonstrate the effectiveness and stability of the level set-based topology optimization framework in handling hyperelastic material distribution.
Following validation of the framework using a cantilever beam benchmark for mean compliance and sensitivity, we now focus on designing a displacement inverter compliant mechanism. This application leverages the large deformation capabilities of hyperelastic materials to minimize mean compliance while achieving the desired displacement inversion.
A 2 × 2m rectangular domain with a thickness of 0.1m is considered for the design. The boundary conditions are illustrated in Figure 6(a). A driving force of f = (0, 12000)N/m is applied at the middle of the left edge, spanning from (0, 0.45)m to (0, 0.55)m. The top section of the left edge (from (0, 0.8)m to (0, 1)m) and the bottom section (from (0, 0)m to (0, 0.2)m) are fully fixed. Additionally, a small resistance force of g = 10N is applied at the middle of the right edge.
Due to the horizontal symmetry of the structure, the design The evolution of the displacement inverter compliant mechanism is illustrated in Figure 7. Starting from an initial uniform material distribution (Figure 7(a)), the topology gradually evolves as material is redistributed to form key structural elements (Figure 7(b, c)). This refinement enhances flexibility and force transmission, ultimately leading to the optimized design in Figure 7(d). The final structure features a lightweight, well-connected topology with smooth boundaries, effectively enabling displacement inversion. The thin central connection functions as a hinge, facilitating large deformations and fully leveraging the hyperelastic material's flexibility, while the oblique beams act as a reinforcing frame, providing structural support with a relatively thicker design.
The convergence history in Figure 8 demonstrates the stability of the optimization process. The objective function decreases sharply in the initial iterations, stabilizes around iteration 36 (Figure 7(c)), and converges at iteration 116. The volume ratio gradually reaches 50%, ensuring material constraint compliance. This demonstrates the efficiency and robustness of the level set-based optimization framework. The finalized 3D structure of the optimized displacement inverter compliant mechanism is shown in Figure 9(b). A prescribed displacement of 0.3m is applied at the middle of the left edge, while the top and bottom sections of the left edge are fixed. The resulting deformation shows that the right end displaces by -0.312m, shown in Figure 9(a), successfully achieving the intended displacement inversion effect. The deformation pattern validates the functionality of the optimized compliant structure, demonstrating the effectiveness of the level set-based topology optimization framework in designing hyperelastic mechanisms
This study presents a level-set-based topology optimization framework for designing compliant mechanisms made of hyperelastic materials, with particular emphasis on neo-Hookean constitutive models. Unlike traditional approaches that focus primarily on stiffness optimization, our framework effectively addresses the challenges of both geometric and material nonlinearities in a unified computational environment. The derivation of shape sensitivity analysis using the adjoint method represents a key advancement, as it successfully incorporates higher-order displacement gradient terms that previous studies often neglected. This comprehensive sensitivity formulation enables more accurate boundary evolution during the optimization process, result-ing in designs that fully exploit the large deformation capabilities of hyperelastic materials. Additionally, our approach eliminates the need for preliminary linear approximations, allowing for direct optimization of nonlinear structures from arbitrary initial configurations.
Through numerical examples, we have demonstrated the effectiveness of the proposed method across different applications. The cantilever beam benchmark validated the accuracy of our sensitivity analysis for mean compliance minimization, showing consistent convergence despite the nonlinear material behavior. More significantly, the displacement inverter compliant mechanism example illustrated how the framework can be applied to design structures that achieve specific motion requirements while leveraging the unique properties of hyperelastic materials.
When compared to topology optimization methods for linear materials, our framework offers several distinct advantages. First, it directly accounts for large deformations in the optimization process, eliminating the geometric limitations typically associated with small-strain assumptions. Second, the level set representation provides clear, smooth boundaries that facilitate manufacturing, particularly important for soft material fabrication where boundary definition is critical. Third, our method naturally handles the complex stress distributions characteristic of hyperelastic deformation, resulting in more reliable performance predictions. Nevertheless, limitations and opportunities for improvement remain. The computational cost of nonlinear analysis within each optimization iteration presents challenges for scaling to very large problems. Future research will focus on extending this framework to address more complex design challenges in soft robotics and biomedical applications. We are developing nonlinear TO methods with three-dimensional volumetric optimization capabilities, incorporating pneumatic actuation mechanisms directly into the formulation, and exploring multi-material optimization for structures with spatially varying stiffness properties. We expect our work can bridge the gap between computational design and physical implementation, enabling a new generation of hyperelastic structures tailored for applications requiring controlled, largescale deformations.
Copyright © 2025 by ASME