Controlling the energy jump of multistable structures using shape optimization

1. Introduction.

• Shape optimization selects a domain to optimize a cost functional subject to physical constraints, often expressed as PDEs.The paper treats domain shape as the control variable in PDE-constrained optimization.Source: §1 · p. 1
• The work targets control of the energy required for metamaterials to snap between stable branches, rather than tuning stress-strain response or negative Poisson’s ratio.Branch switching can dissipate impact energy, making control of that dissipation desirable for impact protection.Source: §1 · p. 1
• Classical shape optimization commonly assumes a unique PDE solution and a single-valued shape-to-solution map, while many physical problems have multiple solutions.For nonconvex energy minimization, these can include local minima, maxima, and saddle points.Source: §1 · pp. 1–2
• Figure 1 illustrates a domain with two holes and three equilibria under the same boundary conditions: two stable states and an unstable state.Dirichlet conditions are imposed on the top and bottom; other boundaries have natural conditions.Source: Figure 1 · p. 2
• The framework parameterizes shape deformation and uses the implicit function theorem to give sufficient conditions for local continuation of solution branches as shape varies.The algorithm tracks those branches throughout shape optimization.Source: §1 · p. 2

2. Problem statement.

• The displacement is a stationary point of potential energy, whose internal energy density is defined using the deformation gradient and neo-Hookean stress.The variational problem is obtained by setting the first derivative of potential energy to zero; P denotes the first Piola–Kirchhoff stress tensor.Source: §2 · pp. 2–3
• The nonlinear variational problem generally admits multiple distinct, isolated solutions on a domain, which are used to define the energy-jump control objective.The problem concerns a neo-Hookean hyperelastic material on a bounded Lipschitz domain in two or three dimensions.Source: §2 · pp. 2–3
• The objective targets signed energy jumps between solution pairs, with pair weights balancing their importance and a penalty preventing distinct solutions from collapsing together.A zero pair weight means that pair’s jump is not controlled.Source: §2 · p. 3
• The objective depends on the domain through integrals and implicitly through its PDE solutions, so the optimization requires suitably regular maps from domains to each tracked solution.Solutions are computed and labelled on the initial domain, then recomputed and labelled after each domain update to maintain consistency.Source: §2 · p. 3

3. Iterative shape optimization.

• The algorithm represents each candidate domain as Ωk = Tk(Ω0) and selects perturbations to produce a decreasing objective sequence.The descent is conditional on suitable selection of perturbations; the paper does not establish that the sequence is minimizing.Source: §3 · p. 4
• The paper focuses on numerical realization and leaves existence and uniqueness for the optimization problem, and conditions guaranteeing a minimizing sequence, unaddressed.The objective sequence is decreasing only under suitable selection of domain perturbations.Source: §3 · p. 4
• Under suitable assumptions, IFT-based continuation consistently pairs indexed solutions across successive domains, supporting objective evaluation for multiple branches.The domain deformation is scaled as a continuation parameter, and sufficiently small updates allow branch correspondence to be checked by solution proximity.Source: §3 · p. 4
• The optimization uses dogleg trust-region steps with BFGS Hessian updates and shape-gradient directions.The implementation uses ROL, Fireshape, Firedrake, pyadjoint, and UFL.Source: §3 · p. 4
• Candidate-shape solutions are computed with Newton iterations initialized from the corresponding solutions on the current domain.When domain updates are sufficiently small, the transferred solution can be a good initial guess and Newton’s method may converge in a few iterations.Source: §3 · pp. 4–5

4. Numerical Results.

• The experiments use perforated elastomeric square domains with N = 2 or 4 hole arrays, and three tracked solutions: two stable branches and an intervening saddle.The domains are also extruded to depth h = 1 for three-dimensional simulations; initial solutions are computed using deflated continuation.Source: §4 · pp. 5–6
• The two objective families target ratios of energy gaps; for the stated symmetric initial branches and central saddle, the adjacent initial gaps are equal.The numerical experiments use Young’s modulus E = 10^6 and Poisson’s ratio ν = 0.3.Source: §4 · pp. 5–6
• For N = 2, Figure 3 shows the initial domain alongside optimized shapes for target energy jump ratios r = 0.6, 0.4, and 0.2, using objectives J1 or J2.J1 appears in the first row and J2 in the second.Source: §4 · p. 6, Figure 3 · p. 6
• For r = 0.4, J2 brings the u1 and u3 energy curves together toward u2, while J1 reduces the u1–u2 gap and leaves u3 uncontrolled.The J1 u3 energy drifts mildly; the J2 pattern is also reported for N = 4.Source: §4 · p. 6, Figure 4 · p. 7
• At r = 0.2, minimizing J2 reaches the target energy-gap ratio, whereas J1 has difficulty and misses it.Figure 5 shows convergence histories for r = 0.2, 0.4, and 0.6 with N = 2 and 4.Source: §4 · p. 6, Figure 5 · p. 7
• For N = 2 and 4, the optimization histories are qualitatively similar, although direct J1 optimization at r = 0.2 fails to reach its target.The authors note that the different J1 and J2 convergence behavior may relate to symmetry preservation or loss, but say further research is needed.Source: §4 · pp. 6–7
• Continuation reaches difficult target ratios, including r = 0.2 and r = 1.6, but requires more iterations across its steps.Direct optimization for large energy-gap increases can be challenging because of competing geometric constraints.Source: §4 · pp. 7–8, Figure 6 · p. 8
• In three dimensions, minimizing J2 reduces the energy gap to approximately 20% of its initial value and stretches the optimized domain opposite the stable-solution branches.The experiment uses N = 2; the optimized geometry is reported to resemble the two-dimensional pattern.Source: §4 · p. 8

• The framework controls energy jumps between distinct solutions of a neo-Hookean metamaterial and may support impact-energy dissipation.Source: §5 · p. 8
• The methodology’s effectiveness is validated through numerical experiments in two and three dimensions.Source: §5 · p. 8
• The method assumes a constant number of solution branches and does not directly apply to self-contact problems.Self-contact replaces the PDE constraint with a quasi-variational inequality.Source: §5 · p. 8

• The deformation gradient is F = ∇φ, with J = det F > 0, and I1 = tr(F^⊤F).Source: Appendix A · p. 11
• The derivation splits the neo-Hookean strain-energy density W into three terms and computes each derivative separately.Source: Appendix A · p. 11
• The contribution of W1 to P is µF.Source: Appendix A · p. 11
• For W3, the directional derivative is λ ln J F^-⊤:H, and P represents the derivative of W at F.Source: Appendix A · pp. 11–12